Computationally efficient and robust kinematic calibration methodologies and their application to industrial robots

Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

Contents lists available at ScienceDirect

Robotics and Computer-Integrated Manufacturing journal homepage: www.elsevier.com/locate/rcim

Computationally efficient and robust kinematic calibration methodologies and their application to industrial robots Temesguen Messay a,n, Raúl Ordóñez a, Eric Marcil b a b

Department of Electrical and Computer Engineering, University of Dayton, 300 College Park, Dayton, OH 45469-0232, United States Yaskawa America, Inc., Motoman Robotics Division, Headquarters and Manufacturing Facility, 100 Automation Way, Miamisburg, OH 45342, United States

art ic l e i nf o

a b s t r a c t

Article history: Received 14 December 2014 Received in revised form 24 May 2015 Accepted 17 June 2015

In this paper, we present new computationally efficient and robust kinematic calibration algorithms for industrial robots that make use of partial measurements. These include a calibration method that requires the supply of Cartesian coordinates of the calibration points (3DCAL) and another calibration technique that only requires the radial measurements from the calibration points to some reference (1DCAL). Neither method requires orientation measurements nor the explicit knowledge of the whereabout of a reference frame. Contrary to most other similar works, both methods make use of a simplified version of the original Denavit–Hartenberg (DH) kinematic model. The simplified DH(-) model has not only proven to be robust and effective in calibrating industrial manipulators but it is also favored from a computational efficiency viewpoint since it consists of comparatively fewer error parameters. We present an analytical approach to develop a set of guidelines that need to be considered in order to properly construct the DH(-) model such that it is parameterically continuous and non-redundant. We also propose an automated method to provide a characterization of the error parameters that is insightful so as to correctly deduce the DH(-) error model of a manipulator. The method makes use of a novel hybrid optimization scheme to conduct a statistical analysis of the error parameters that is indicative of their relevance. We made note that, for the industrial robots used in this paper and similar ones, calibrating the home position only is sufficient to attain adequate results for most robotic applications. Hence, we put forward for consideration a yet simpler calibration model; the DH(-)(-) model. We employ the Trust Region (TR) method to minimize the objective functions of both frameworks (3DCAL and 1DCAL). The performance of the proposed methods is compared to that of a state-of-the-art commercial system (MotoCal) using the same materials, data and internationally recognized performance standards. Our experimental results suggest that our methods yield improved results compared to that of MotoCal. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Kinematic calibration Optimization Simulated annealing Trust region CompuGauge MotoCal Industrial robots Yaskawa Motoman Robotics, Inc

1. Introduction Repeatability (also known as test–retest reliability) and accuracy are important characteristics of industrial robots [1–5]. The demand of improving those two qualities of industrial robots has therefore been growing continuously over the past two decades. It is a known fact that today's industrial manipulators have satisfactory repeatability (better built) but poor accuracy due to numerous sources of errors [3,6–9]. Robot calibration is the process of enhancing the positioning accuracy of a manipulator through software rather than modifying the mechanical structure (i.e., design) of the robot itself [10–12]. The leading source of lack of accuracy is the discrepancy between the mathematical model of n

Corresponding author. E-mail addresses: [email protected] (T. Messay), [email protected] (R. Ordóñez), [email protected] (E. Marcil). http://dx.doi.org/10.1016/j.rcim.2015.06.003 0736-5845/& 2015 Elsevier Ltd. All rights reserved.

the manipulator in the controller and the actual geometry of the structure [2,3,6,8,13,14]. An accurate representation of the geometry of a robot is very crucial because the efficiency of methods for planning and controlling robot motions is highly dependent on such a mathematical model [12,15–19]. Hence, in this paper, we present computationally efficient and robust kinematic calibration techniques for industrial robots. A kinematic calibration procedure involves modelling, measuring, identifying parameters, implementing compensation and validating [10,20,21]. 1.1. Modelling The first procedural step is to derive a mathematical model that relates the robot's joint angles to the pose of its end-effector and that takes into account the geometric error parameters that need to be identified. The standard Denavit–Hartenberg (DH) convention [22] is universally used for kinematic modelling in robotics [3,22–24,16,18]. However, a standard DH based kinematic

34

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

calibration model is singular for manipulators with two consecutive parallel (or near parallel) joint axes [25–27]. Since the majority of industrial robots possess at least two parallel joint axes, significant efforts were made to solve such a problem. Authors either proposed to use a modified version of the standard DH convention (e.g., adding an extra parameter to the original DH model) or introduced their own model to resolve the presented challenge [10,20,25,26,28–33]. However, a potential drawback of some of those alternate solutions is that the compensation is not directly implementable in controllers of existing industrial manipulators. In this paper, we propose to “simplify” the standard DH based kinematic calibration model so as to make it continuous and nonredundant (i.e., we opt to not consider some of the DH error parameters deemed redundant and hence irrelevant). Our approach has not only proven to overcome the problem of discontinuity [3,6,25–27] and redundancy [3,6] but it is also favored from a computational efficiency viewpoint since it consists of comparatively fewer error parameters. Furthermore, we believe that it is applicable and well-suited to calibrate industrial manipulators like those of Yaskawa Motoman Robotics, Inc. (or other similar manipulators) because not all DH parameters are accessible to be modified [34–36]. We have established useful rules/ guidelines based on empirical studies and generalized geometric knowledge of kinematic modelling that need to be considered in order to make the error model non-redundant. We shall show that the resultant simplified DH(-) model is more robust compared to existing calibration methods as it ensures that a reliable end-result is attained. We have also devised an automated method to provide a parameter assessment that is greatly insightful in identifying irrelevant parameters (i.e., that is helpful to properly construct the DH(-) model of a given manipulator). In particular, we believe that this automated assessment to be greatly advantageous when calibrating manipulators of complex geometry (e.g., structures with higher Degree of Freedom (DOF) compared to the ones considered in this work). The method makes use of a novel hybrid optimization scheme composed of the Simulated Annealing (SA) algorithm [37–41] and the Trust Region (TR) algorithm [42–46]. It provides a statistical analysis on the estimates of a given error parameter that is suggestive of its relevance. Thus, for the error parameters that the auto-rating is above a given threshold or the end-user so designates, only a basic understanding of the geometry of the robot is required from the user to determine which ones are nonpertinent (i.e., that consequently need to be discarded in order to correctly derive the DH(-) kinematic calibration model for the manipulator at hand). Additionally, we shall demonstrate that, for the type of robots (e.g., [47–49]) and controllers (e.g., [34–36]) that are put to use in this work, precision inaccuracy is mostly due to incorrectness of the offset values used to describe the manipulator's home position (also known as Absolute [ABS] data [34–36]). Hence, we further simplify the DH(-) model and introduce the DH(-)(-) model, a simpler model that consists of calibrating the home position (and tool if applicable) only. We shall demonstrate that, although incomplete, the DH(-) and DH(-)(-) kinematic calibration models are capable of exceptional performance and can be used to calibrate a wide range of industrial robots provided that the models are accurately constructed. 1.2. Measuring In practice, another crucial choice is the measurement system. The efficiency of the calibration method is a function of the accuracy of the observations (measurements) itself [3,6,50,51]. A number of different measurement systems have been used for

robot calibration and/or validation. Generally speaking, the measurement systems used to calibrate a robot can be classified into two groups: complete pose measurement and partial pose measurement. A complete pose measurement of the tool pose would consist of three position coordinates and three orientation angles (6D). This type of measurement method yields the maximum information for a given configuration. Examples of kinematic calibration techniques that make use of complete pose measurement include those described in [6,10,51–57]. Partial measurements of the tool pose (that is, less than six measured values per observation ranging from 3D to 1D) can also be used to identify the robot's geometry. Methodologies that have employed partial measurement equipments to calibrate a robot can be found in [3,20,23,24,26,50,58–74]. What is paramount to note from the cited research efforts (type of measurements used in their methods and reported performance analyses) is that incomplete pose information in accordance with a specific data-collection scheme can be used to successfully calibrate a robot (yields satisfactory results for most industrial robotic applications). Also to note is that, in most cases, complete pose measuring devices are found to be relatively more expensive. For these reasons, it behooves us to develop a kinematic parameters identification solution that functions appropriately with incomplete pose information. Moreover, considering the variability among existing partial pose measuring devices, we believe that, in order for a solution to have practical relevance, it ought to be adaptable to the supply of various types of incomplete pose information. Hence, our overall kinematic calibration framework is composed of two methodologies: one that uses 3D position measurements (3DCAL) and another that requires radial distance measurements only (1DCAL). Note that neither method requires orientation measurements. In this work, we made use of the 3D and 1D CompuGauge [75] measurement systems to acquire the partial pose measurements of the calibration points. 1.3. Identifying parameters For the third procedural calibration step, in our framework, we have formulated kinematic parameters identification process as a simplified optimization problem. In our 3DCAL and 1DCAL systems, we make use of the DH(-) and DH(-)(-) error models. Unlike the 3D calibration methodologies presented in [3,58] and the 1D calibration scheme proposed by [71], our approaches do not require the explicit knowledge of the position of a reference frame (neither the robot's base frame nor the whereabouts of the measurement equipment's reference frame). We therefore believe that our solutions are more flexible compared to that of [3,58,71]. In this study, the TR approach and the Levenberg–Marquardt (LM) algorithm [76–78] are implemented to minimize the objective functions of both systems (3DCAL and 1DCAL). Even though both (TR and LM) algorithms converged to the same numerical values for all cases of non-linear least square analyses, we found TR to be more efficient compared to LM and hence the former is favored and used to optimize over the cost functions of both systems (i.e., used to solve for the relevant error parameters of the DH(-) and DH(-)(-) error models of both systems). 1.4. Implementing compensation For this calibration procedural step, we allow user interaction. It must be recalled that access to modify the parameters in some controllers is not obvious, and may not be possible for second and third parties. Hence, we allow the user to select which DH parameters estimate subject to our rules that need to be applied to successfully derive the DH(-) error model (i.e., in order to render the original DH convention based error model continuous and

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

non-redundant). It is important for the compensation of a kinematic calibration technique to be practical that it is transferable to control software of existing industrial manipulators [10]. From that perspective (i.e., considering the widespread use of the original DH convention in existing robot control software), we believe that our approach is desirable for this procedure since DH parameters are obtained directly from both methods (3DCAL and 1DCAL). 1.5. Validating: performance analysis and comparison It is a difficult task to make definitive comparison against and/ or among previously published calibration systems due to variability in materials, error models and performance criteria used in the methodologies. It is well known that the performance results of any calibration system can differ significantly depending on those variables. Fortunately, a comprehensive comparison is made possible here thanks to Yaskawa Motoman Robotics, Inc., who allowed us to make use of CompuGauge [75] and MotoCal [79] (commercial systems that are currently used in their setting to calibrate and authenticate various manipulators). Our experimental results are compared against those of MotoCal [79]. We have used the same material and performance metric to put our performance into clear context. We shall show that our calibration methods are efficient and yield improved results compared to the state-of-the-art system of the same genre. 1.6. Paper organization The remainder of this paper is organized as follows. We begin by describing the materials that we used to carry out numerous experiments in Section 2. In Section 3, we present two reliable techniques that can be used to properly construct the DH(-) and DH(-)(-) error models. The 3DCAL and 1DCAL calibration schemes are also described in that section. Experimental results and related discussion are presented in Section 4. Finally, in Section 5, we offer conclusions. Where relevant, some of the previously published calibration methods are discussed.

2. Materials 2.1. Manipulators: the HP20D and MH6 To test our 3DCAL system, we made use of the Yaskawa Motoman Robotics Inc. HP20D manipulator [47] and the control software of the DX100 controller [35]. The HP20D is a versatile high speed 6 DOF robot with all six revolute joints [47]. It is an anthropomorphic robot composed of a spherical wrist [16,18,12]. The HP20D is used in a variety of applications such as cutting, dispensing, handling, machine tending, material removal, packaging, press tending and welding [47]. Details concerning the specifications of the robot and controller are presented in [47,35]. The DH Table according to the convention described in [16] is presented in Table 1. Table 1 DH parameters of the HP20D manipulator. Link

d i in mm

a i in mm

α i in degrees

θ i in degrees

1 2

0 0

150 760

90 0

3 4 5

0 795 0

140 0 0

90 90 90

6

105

0

0

θ1 θ 2 + 90 θ3 θ4 θ5 − 90 θ6

35

Table 2 DH parameters of the MH6 manipulator. Link

d i in mm

a i in mm

α i in degrees

θ i in degrees

1 2

0 0

150 614

90 0

3 4 5

0 640 0

155 0 0

90 90 90

6

95

0

0

θ1 θ 2 + 90 θ3 θ4 θ5 − 90 θ6

The MH6 is a similar 6DOF serial robot that is used for alike applications as the HP20D [48]. We made use of the MH6 manipulator and DX100 controller to validate our 1DCAL method [48,35]. Specifications of the MH6 robot can be found in [48,35]. Its DH Table is shown in Table 2. 2.2. Measurement systems: 3D and 1D CompuGauge For the purpose of measuring the Cartesian coordinates of the calibration points for our 3DCAL system, we made use of 3D CompuGauge [75]. CompuGauge is a system that allows 3D measurement and performance testing of robot positioning. It is composed of hardware and software. As shown in Fig. 1, the hardware consists of two triangulation beams. The “measurement attachment” (tool zoomed-in in Fig. 1) is mounted on the robot. The four “measurement cables” originating from both triangulation beams are connected to the tool. High resolution, low inertia optical encoders are used to constantly measure the extension of the cables. The data is interpreted by the software to deduce an accurate measurement of the test point with respect to the user-defined reference frame [75]. Note that only three out of the four cables that originate from the fixed triangulation beams are utilized to calculate the coordinates of a given test point. The extra cable is dedicated for verification purposes. In addition to data acquisition, the software can also perform various robot performance tests following two different internationally recognized performance standards: “Manipulating Industrial Robots – Performance Criteria and Related Testing Methods” (ISO 9283) and “American National Standard for the Evaluation of Point-to-Point and Static Performance Characteristics of Industrial Robots and Robot Systems” (ANSI/RIA R15.05) [75]. In this paper, we make use of that software capability to measure and report our performance. More will be said about the performance metric and analysis in Section 4. Specifications of the discussed measurement system are presented in Table 3. In order to acquire the 1D radial measurements of the calibration points for our 1DCAL system, we make use of 1D CompuGauge. 1D CompuGauge consists of an encoder, a specially designed end-effector that is to be mounted on the robot and software. In essence, it is a simplified version of 3D CompuGauge. As shown in Fig. 2, the encoder is to be anchored outside of the manipulator's workspace; preferably not greater than 45 cm to any of the two sides of the robot's base and approximately the same height vis a vis to the base. This location permits the largest working envelope, reduces interference with robot components and increases calibration accuracy. The software provides similar capabilities as the ones of 3D CompuGauge; a means for data acquisition and robot performance testing based on alike standards. 3. Proposed calibration methods: 3DCAL and 1DCAL In this section, we describe the proposed calibration methods. The analytical approach used to develop a set of guidelines that

36

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

Fig. 1. The 3D CompuGauge measurement system and attachment (tool) [75]. Table 3 Specifications of the 3D CompuGauge system [75]. Performance Resolution Repeatability Accuracy Measurement space Tracking rate Sampling frequency Data range Measurement Time

0.010 mm 0.020 mm 7 0.150 mm 1500 mm  1500 mm  1500 mm up to 5 m/s from 25 to 1000 Hz up to 32000 points/file up to 10 min

Mechanical (per “triangulation beam”) Dimensions (approx.) Weight (approx.):

(L ) 1730 mm × (W ) 165 mm × (H ) 165 mm 7 kg

Environmental Operating temperature:

20 7 2 °C (687 5)

need to be applied in order to adequately construct the DH(-) and DH(-)(-) kinematic calibration models are discussed in Section 3.1. 3DCAL is described in Section 3.2. 1DCAL is presented in Section 3.3. In Section 3.4, two optimization techniques for the purpose of minimizing the cost functions of 3DCAL and 1DCAL (identifying the error parameters of the non-redundant DH(-) and DH(-)(-) kinematic calibration models) are compared. Following that, an automated parametric relevance assessment technique that can greatly aid the end-user in deducing the DH(-) model of a given manipulator is presented in Section 3.5. In that section, we also correlate the results of the analytical and automated approaches in order to demonstrate the agreement among the two approaches in finding the simplified DH(-) model. 3.1. The DH(-) and DH(-)(-) kinematic calibration models: the analytical approach A kinematic calibration model is a mathematical function that relates the robot's joint angles to the pose of its end-effector and

Fig. 2. 1D CompuGauge System [75,79].

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

that takes into account the geometric error parameters that need to be modeled [12,15–18]. As told in Section 1, the standard DH model is currently and widely used in robotics and hence is our first candidate under consideration. Hayati et al. [25,26] first made note that a standard DH model based kinematic calibration scheme is not parameterically continuous. More specifically a small misalignment of parallel joint axes, which corresponds to error in nominally zero valued DH parameter α (also known as link twist [12,15–18]), causes a large change in the other DH parameters [25,26]. Notwithstanding that Nubiola and Bonev [3,6] recently proposed an elegant solution to overcome the problem of disproportion at parallel or near parallel joint axes. In their method, they forced axis 3 to be parallel to axis 2 for simplicity (i.e., they chose not to perturb the zero valued DH parameter α) and were able to successfully calibrate a manipulator by making use of the classic DH model. In this work, we apply that solution in order to meet the “Continuity” criterion described in [27,32] and make use of an original DH convention based kinematic calibration model to calibrate various manipulators. Table 4 presents the continuous kinematic calibration model applied to the HP20D manipulator. The HP20D robot's kinematic calibration model can be mathematically represented as

T (q, e)0tool = A 1(q1, e1) A 2 (q2, e 2 ) A 3(q3, e 3 ) A 4 (q4 , e 4 ) A 5, e5 (q5 ) A 6 (q6, e 6 ) A tool (e tool ),

(1)

where q is a vector containing the joint values (i.e., [q1, q2, … , q6 ] which in our case = [θ1, θ 2, … , θ 6 ] since we are dealing with 6DOF manipulator that possesses revolute joint axes). Vector e is composed of the sub-vectors [e1, e2, … , e6 , e tool ]. More specifically, the sub-vector e i contains the deviations δdi , δai , δαi and δθ i that are to be added to the corresponding nominal values where applicable (i.e., as shown in Table 4). Similarly, sub-vector e tool has the error parameters δXtool , δYtool, δZtool , δRXtool, δRYtool, δRZtool as contents. It follows that, A i (qi , e i ) and A tool (e tool ) can be calculated using Eqs. (2) and (3), respectively. That is,

A i (qi , e i ) = Rot (z, θ i + δθ i ) Trans (z, di + δdi ) Trans (x, ai + δai ) Rot (x, αi + δαi ),

(2)

A tool (e tool ) = Trans (x, Xtool + δXtool ) Trans (y , Ytool + δYtool ) Trans (z, Ztool + δZtool ) R tool,

(3)

where

R tool = RPY (RXtool + δRXtool, RYtool + δRYtool, RZtool + δRZtool ).

(4)

transformation that takes into account the geometric error parameters; A i (qi , e i )). Note that 6 out of the 28 error parameters of the presented DH based calibration model are dedicated for estimating the position and orientation of the end-effector itself. Let p = [x, y, z ]T denote the calculated position of the end-effector. In other words, the vector p contains the first three elements of the last column of T (q, e)0tool for some configuration q and a specified set of error parameters e. The robot's DH kinematic calibration model can be expressed in a compact form as

p = f (q, e),

(5)

where vector e belonging to Ω contains all 28 geometric error parameters and Ω is the set of allowable deviations (i.e., for a given error parameter e, e ∈ [ − e, e¯ ] where − e and e¯ are robot manufacturer dependent scalar range limits). More will be said on how to define the space Ω in Section 4.3. We believe that forcing the alignment of parallel or near parallel joint axes as done by Nubiola and Bonev [6,3] is well applicable in our case because access to modify the DH parameter α in Yaskawa Motoman Robotics, Inc. controllers is not possible [34– 36]. Hence, in addition to that, we extend the idea and elect to “freeze” all DH α parameters (i.e., making use of the nominal values) in our calibration schemes. Moreover, by examining the rank of the Jacobian of their leastsquare optimization method, Nubiola and Bonev [3,6] found that several error parameters are often redundant. In order to attain a stable solution, they chose to take out what they referred to as “useless” error parameters [6,3]. For example, the error parameters δd2 and δd3 are dependent because of the parallel axes of joints 2 and 3 (which are forced to be perfectly parallel) and so they chose to not consider δd3 in their model [3,6]. Another form of redundancy that can be problematic in achieving convergence (stable solution) for an optimization scheme like the one we will be employing can be demonstrated using the example shown in Fig. 3. The figure shows that different combinations of the other three parameters θi, ai, and di lead to the same transformation from point A to B. In order to circumvent the discussed issues in the previous paragraph, our product usage guidelines include the following terms. Since axes 2 and 3 of the HP20D are considered perfectly parallel, only one of the two DH parameters that describes an offset along that dimension (e.g., δd2 or δd3) is to be considered in the model. Moreover, for the reason explained in the above, one of the two DH parameters (ai or di) of a transformation that describes the spatial relation of any two given links (i.e., A i ) must be kept “frozen” at its nominal value. In fact, in this work, we will not be considering zero valued DH parameters (ai and di) in our kinematic

RPY in Eq. (4) is a function that, given the Roll/RZtool, Pitch/RYtool and Yaw/RXtool angles of the end-effector mounted on the manipulator, returns a 4  4 homogenous transformation matrix R tool representative of the orientation of the tool frame [15]. See Appendix A for details of the abbreviated homogeneous transformation matrices (including the explicit expression of a spatial Table 4 Original DH convention based continuous kinematic calibration model applied to HP20D manipulator. Link

d i in mm

a i in mm

α i in degrees

1

0 + δd1

150 + δa1

2

0 + δd2

760 + δa2

90 + δα1 0

3

0 + δd 3

140 + δa3

90 + δα3

4

795 + δd4

0 + δa 4

−90 + δα4

θ 4 + δθ 4

5

0 + δd5

0 + δa 5

(θ5 − 90) + δθ5

6

105 + δd6

0 + δa 6

90 + δα5 0

37

θ i in degrees θ1 + δθ1 (θ 2 + 90) + δθ 2 θ 3 + δθ 3

θ 6 + δθ 6 Fig. 3. Illustration of DH parameters redundancy.

38

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

calibration models. The reason for that lies in a long standing experience (achieved over the course of calibrating numerous Yaskawa Motoman Robotics, Inc. manipulators) and extensive experiment done by us, that including those particular parameters does not contribute a noteworthy overall performance gain in return. Moreover, note that this general rule also complies with the guidelines established to prevent redundancy (i.e., only one of the two ai or di can be calibrated) because, as it can be seen in Tables 1 and 2, the di parameters of the arm of those robots are all zero valued while the ai parameters are nonzero elements. Likewise, all of the ai parameters related to the spherical wrist of the HP20D and MH6 are zero valued and most of the corresponding di parameters are nonzero in value. We believe that this general rule (of not calibrating zero valued parameters) can also be extended to other industrial robots that exhibit similar anthropomorphic design (e.g., the Unimation PUMA 560). However, if the design of the manipulator greatly differs compared to the ones put to use in this work (or similar ones), a revision to this general rule may be in order. Furthermore, we will not be considering δθ1. As mentioned in Section 1, our approach is invariant to the fixed frames (robot's and measurement equipment's reference frames) as the Euclidean distances among the calibration points are the measurements put to use. Hence, estimating that particular error parameter will only translate and/or rotate all calibration points while the Euclidean distances among those observations remain unchanged. Since that will bring about undesired redundancy, in our approach we shall not perturb the error parameter associated with the first joint angle θ1 during calibration as well. Finally, another aspect of kinematic calibration is the need to determine the tool transform when necessary. The provided measurement-attachment/tool shown in Fig. 1 is well defined and ideally is to be mounted directly onto the robot. Unfortunately, for the HP20D, we used a fabricated mounting plate due to hardware incompatibility and hence the error parameters of A tool (e tool ) in Eq. (3) need to be identified. More will be said about the experimental procedure that we performed to obtain an initial tool definition in Section 4. What needs to be emphasized here is the guidelines that need to be set in order to attain a good estimate after an initial guess of the tool specifications. It follows that if the attached tool is well defined, then error parameters listed in the last/sixth row of Table 4 (δd6 and δθ 6 ) can be estimated. However, if one needs to estimate tool definition (as done in this work for the case of the HP20D robot) then those error parameters must not be taken into account for the obvious reason to avoid redundancy. This will be justified via our auto-assessment technique in Section 3.5 and further discussed in Section 4. In the case of the HP20D, we made use of a flat fabricated mounting plate (approximately 20 mm in thickness), one that guarantees certainty of the orientation of the tool frame and hence we decided to estimate the error parameters related to the position of the tool frame only (δXtool , δYtool and δZtool ). Let “DH(-) model” denote the kinematic calibration model after modifying the model described in Table 4 as discussed thus far (i.e., after applying the discussed guidelines). Table 5 presents the DH(-) kinematic calibration model for the HP20D manipulator that we will use to evaluate our 3DCAL methodology. The “simplified” model is now composed of only 11 error parameters (where 3 parameters are dedicated for estimating the location of the tool frame). Note that although our kinematic calibration model is incomplete (i.e., does not meet the “completeness” criteria defined in [27,80]), we shall show that it is capable of exceptional performance (i.e., yields accuracy that is more than satisfactory for most robotics applications). More will be said about this in Section 4. Also to note is that the DH(-) model is without a doubt favorable (compared to that of the original DH or other complex models

Table 5 DH(-) kinematic calibration model of the HP20D manipulator. Link

d i in mm

a i in mm

α i in degrees

1

0

150 + δa1

90

θ1

2

0

760 + δa2

0

3

0

90

4

795 + δd4 0 105

140 + δa3 0

(θ 2 + 90) + δθ 2 θ 3 + δθ 3

0

90

0

0

5 6

90

θ i in degrees

θ 4 + δθ 4

(θ5 − 90) + δθ5 θ6

Table 6 DH(-) kinematic calibration model for MH6 manipulator. Link

d i in mm

a i in mm

α i in degrees

1

0

150 + δa1

90

θ1

2

0

614 + δa2

0

3

0

90

4 5

640 + δd4 0

155 + δa3 0

(θ 2 + 90) + δθ 2 θ 3 + δθ 3

0

90

6

95 + δd6

0

0

90

θ i in degrees

θ 4 + δθ 4

(θ5 − 90) + δθ5 θ 6 + δθ 6

discussed in Section 1) from a computational efficiency stand point of view since it consists of fewer degrees of error parameters. In order to validate our 1DCAL approach we made use of the DH(-) model of the MH6 manipulator and the 1D CompuGauge measurement equipment. The continuous and non-redundant DH (-) model for the MH6 manipulator can be constructed by applying the same guidelines as for the HP20D and is presented in Table 6. We have decided to calibrate θ6 and d6 in the case of the MH6 because we made use of a well defined tool that we directly mounted on the robot (i.e., we do not take into account the error parameters related to the tool data). A tool can be accurately calculated using

A tool = Trans (x, 0.951) Trans (y , − 149.920) Trans (z, 138.040) R tool,

(6)

where R tool is a 3  3 identity matrix. All dimensions in that equation are metric (mm). Note that the DH(-) model of the MH6 consists of only 10 error parameters. As briefly mentioned in Section 1 and as it will be demonstrated in Section 4, when it comes to industrial manipulators like the ones of Yaskawa Motoman Robotics, Inc. and alike a great amount of positioning inaccuracy is mostly due to erroneous offset values used to describe the manipulator's home position (also known as ABS data in controllers found in [34–36]). We therefore further simplify the DH(-) model and introduce the DH(-)(-) model, where error parameters related to the home position and tool (if applicable) are only estimated. Based on the logic deliberated in the passage above concerning the error parameter δθ1 (i.e., since our calibration framework is weakly dependent on that parameter), it will not be taken into account in this simpler error model as well. We acknowledge that the newly introduced model is highly incomplete (according to [27,80] and compared to the DH (-) model) but it produces results that we believe are valuable for a majority of robotic applications. Table 7 presents the DH(-)(-) kinematic calibration models for the two manipulators (the HP20D and the MH6, respectively). Note that the DH(-)(-) error model of the HP20D is a function of 7 error parameters (where 3 are tool frame related estimates since we chose to calibrate the tool of vaguely known dimensions). In contrast, the DH(-)(-) model of the MH6 is a function of only 5 error parameters since we chose not to calibrate the tool definition. Also to note is that DH(-)(-) model

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

CompuGauge measurement equipment (i.e., a set of measured Cartesian coordinate vectors and the distinct non-zero values of their Euclidean distance matrix, respectively). As shown in Fig. 4, the multi-variable objective function of our 3DCAL system that is to be minimized (off-line) using an optimization scheme can be expressed as

Table 7 DH(-)(-) kinematic calibration models. Link

d i in mm

a i in mm

α i in degrees

θ i in degrees

The DH(-)(-) kinematic calibration model of the HP20D Manipulator 1 0 150 90 θ1 2 0 760 0 (θ 2 + 90) + δθ 2 3 0 140 90 θ 3 + δθ 3 4 795 0 90 θ 4 + δθ 4 5 0 0 90 (θ5 − 90) + δθ5 6

105

0

0

y3DCAL (e) =

95

0

0

θ6

(N ×(N − 1)) /2

∑ i=1

e⁎ = arg min y3DCAL (e). e∈ Ω

θ 6 + δθ 6

with notable fewer parameters can be computationally processed a lot quicker compared to any other kinematic calibration models. 3.2. 3DCAL This subsection describes the 3DCAL technique that is employed to identify the error parameters of the DH(-) and DH(-)(-) models of the HP20D manipulator. Fig. 4 shows a block diagram of the calibration algorithm where the HP20D manipulator and its DH(-) error model are used as an example. A natural extension of Eq. (5) described in Section 3.1 is

P = f (Q, e),

1 (N × (N − 1))/2)

Δli2,

(8)

∼ where Δli = li − li represents the residual between the calculated and measured distances of the ith Euclidean distance (for a total of (N × (N − 1)) /2 unique Euclidean distances). The optimal perturbation vector e⁎ can then be calculated as

The DH(-)(-) kinematic calibration model of the MH6 Manipulator 1 0 150 90 θ1 2 0 614 0 (θ 2 + 90) + δθ 2 3 0 155 90 θ 3 + δθ 3 4 640 0 90 θ 4 + δθ 4 5 0 0 90 (θ5 − 90) + δθ5 6

39

(7)

where the matrix P is now composed of a set of calculated Cartesian coordinate vectors ⎡⎣p1(q1, e) , p2 (q 2, e) , … , pN (q N , e) ⎤⎦ of N calibration points, and matrix Q contains the corresponding set of joint values (i.e., set of N configuration vectors). The Euclidean distance matrix D of matrix P is calculated [81– 84]. Matrix D is a hollow (diagonal elements are all equal to zero) and symmetric matrix [81–84] and is of size N  N for a given set of N calibration data points. Let l denote a vector that contains the (N × (N − 1)) /2 unique (distinct) non-zero values of matrix D (distinct Euclidean length measurements). Furthermore, let matrix ∼ P and vector ˜l represent arrays that contain the corresponding measured values of those calibration points acquired using the 3D

(9)

Note that the error function y3DCAL (e) that is to be minimized is based on the discrepancy between the calculated and measured Euclidian distances among the tool end-points for the designated configurations. In other words, we do not explicitly make use of the Cartesian coordinates of the calibration points with respect to a fixed frame (neither the robot's base nor 3DCompuGauge's reference frames). We believe that our approach is therefore more practical and robust compared to the ones described in [3,58]. Those two methods call for an extra procedure to determine a common reference coordinate system beforehand and that extra step alone could introduce errors that could potentially propagate down to the end-results [10]. At this point, the challenge is to optimize over the cost function described in Eq. (8) (i.e., calculating the optimal perturbation vector e⁎ as denoted in Eq. (9)). It is basically a nonlinear leastsquares optimization problem that can be solved using either the Trust Region (TR) algorithm [42–46], or Levenberg–Marquardt (LM) algorithm [76–78]. The two are powerful optimization algorithms for large-scale nonlinear problems. In Section 3.4, we will compare the two and make a recommendation with regard to which one is more effective and efficient in identifying the optimal perturbation parameters of the DH(-) and DH(-)(-) error models of the HP20D manipulator using the discussed 3DCAL method. 3.3. 1DCAL The system presented in this subsection requires the supply of radial measurements of the calibration points with respect to

Fig. 4. 3DCAL block diagram. The DH(-) model of the HP20D is used as an example.

40

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

Fig. 5. 1DCAL block diagram. In this figure, the DH(-) model of the MH6 is used as an example

some fixed reference. Fig. 5 presents a block diagram of the 1DCAL approach where the MH6 manipulator and its DH(-) error model are used as an example for illustration purposes. Our 1DCAL system is a simplified version of 3DCAL and builds on the work of Goswami et al. [71]. However, we depart from that in [71] and we initially estimate the location of the fixed reference and then iteratively refine that finding along with the selected DH error parameters via our algorithm during the calibration process. As illustrated in Fig. 5, we propose to define the main objective function that is to be minimized to estimate the location of the fixed reference (encoder) pe = [x e , ye , ze ]T , and the optimal vector e that contains the selected DH error parameters as

y1DCAL (pe, e ) =

1 N

N

∑ (ri − ∼ri)2, i=1

(10)

where ri and ∼ ri are the computed and measured radial distances at the ith calibration point, respectively (for a total of N calibration points). More specifically, ri can be calculated as

ri = ∥ pe − p i ∥,

(11)

where the ith vector pi that describes the position of the end-effector can in turn be calculated by invoking Eq. (5) (i.e., pi = f (q i, e )). This is another nonlinear least-squares optimization problem where either of the two powerful optimization algorithms mentioned in Section 3.2 can be used to optimize over Eq. (10) such that a minimizer

e⁎ = arg

min y1DCAL (pe, e ),

e ∈ Ω, p e ∈ V

(12)

is found. The bounded space V in Eq. (12) is a spherical volume. More will be said about the volumetric space V and the set of allowable deviations Ω in Section 4.3. Note that our approach consists of solving for both minimizers e⁎ and p⁎e simultaneously (although the former vector contains the relevant information to calibrate the manipulator). Also to note is that those two vectors belong to two drastically different subspaces (i.e., Ω and V, respectively). We believe a zero valued vector ∈Ω is an appropriate initial-condition for the purpose of solving for e⁎ since it contains deviations of small magnitude. However, the vector pe contains the location of the encoder and hence a zero

valued vector ∈V is far from the minimizer p⁎e . Thus, to facilitate computation, we take on the task of supplying an adequate starting-point in order to solve for p⁎e more efficiently. We initially optimize over

y1DCAL (p e, 0 ) =

1 N

N

∑ (ri − ∼ri)2,

(13)

i=1

so as to find the minimizer

p⁎e0 = arg

min

0 ∈ Ω, p e ∈ V

y1DCAL (p e, 0 ),

(14) p⁎e .

that is relatively closer to the solution Hence, as an initial procedure, we execute Eq. (14). We then make use of the obtained result for p⁎e0 as an initial starting point to optimize over Eq. (10) so as to solve for e⁎ and p⁎e simultaneously. Once more we are in need of an effective optimization scheme to estimate the error parameters of the MH6 using 1DCAL. In the next subsection, we compare the TR and LM optimization schemes and provide our suggestion concerning which of the two is preferable in terms of effectiveness and efficiency to solve the non-linear least squares optimization problems that we have formulated thus far. 3.4. Optimization scheme: identifying error parameters In this subsection, we discuss the optimization scheme employed to solve for the optimal perturbation vector of the simplified models (DH(-) and DH(-)(-)) using the 3DCAL and 1DCAL frameworks. Aforementioned, both Eqs. (8) and (10) are nonlinear least-squares optimization problems that can be solved using either Trust Region (TR) approach [42–46] or Levenberg–Marquardt (LM) method [76–78]. TR and LM are two efficient optimization algorithms that can handle large-scale nonlinear optimization problems. The details of the two optimization algorithms can be found in [42–46,76–78,85]. LM and TR are both restricted Newton step based methods and hence exhibit quadratical speed of convergence for iterations near the solution e⁎. Moreover, they both can handle singular system of non-linear equations. However when we are far from the solution, LM will slow down drastically since it relies mostly on the gradient search method. In contrast, the solution trajectory of TR can

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

potentially perform a very long leap and expedite the search process. This efficient facet of TR has been proven by numerous literatures including [85,86]. In particular, as stated by [85], “[TR] methods are an evolution of [LM] algorithms. [TR] methods are able to follow the negative curvature of the objective function. [LM] algorithms are NOT able to do so and are thus slower.” We strongly concur with these conclusions as TR demonstrated a slightly faster convergence rate for all cases of nonlinear least square minimization problems presented in this paper. These include Eqs. (8), (10) and (13). For example, in the case of minimizing Eq. (13), given the same initial condition (a zero valued pe0 vector) and identical tolerance based stopping-criterion, LM carried out 13 iterations while TR only performed 7 iterations to converge to the same solution. Hence, even though they both converged to the same numerical solutions for all calibration cases, we favor TR (in lieu of LM) as it is computationally more efficient for this application. In this work, the termination condition for TR that we used is that the deviation error is small enough. To be precise, until the deflection error is in the order of 10 12. We also made use of TR as a solver for the purpose of carrying out a parametric redundancy auto-assessment using a Simulated Annealing (SA) based approach that we are about to introduce in the next subsection. 3.5. The DH(-) kinematic calibration model: the automated method Here, we propose an automated method that calls for numerical analysis to provide a redundancy assessment of the error parameters. As mentioned in the previous subsection, TR is found to be adequate and efficient to identify the error parameters of the DH(-) and DH(-)(-) kinematic calibration models (i.e., solving nonredundant systems). However, when it comes to redundant error models (e.g., original and continuous DH based error model of the HP20D manipulator shown in Table 4), TR and similar solvers are sensitive to the supplied initial condition when searching for the optimizer(s) as error models with redundancy are problems that contain multiple minima. We propose to take advantage of such sensitivity in order to identify redundant parameters by running multiple searches. We shall examine multiple initial starting points using TR and inspect if the selected error parameters converge to the same numerical values in which case those will be deemed non-redundant. Conversely, error parameters that converge to different numerical values but without any significant impact in the output of the objective function are adjudged ineffectual, and hence, redundant. We propose to make use of another type of solver to gather a rich set of starting points. We then run TR using those starting points in order to perform a search for multiple optimizers. In essence, we have a global search mechanism similar to the ones described in [87–89] with the distinct exception that we make use of SA to generate starting point candidates. SA is a stochastic optimization algorithm on its own right [37– 41]. Note that in this work it was considered to use SA as a stand alone solver for this application (to identify the error parameters of the DH(-) and DH(-)(-) models) but we found SA to be inadequate in converging to a useful result. Difficult and arduous tasks that need to be mastered for SA to be effective include a proper choice of initial temperature and temperature model for each dimension of the solution space, setting tuning parameters (e.g., annealing constant and re-annealing interval), determining search stopping criteria (e.g., maximum number of iterations, change in function or parameters value [tolerance], and time allowed) [37–41]. However, we made note of SA's potent capability of intelligently exploring the solution space in question. It is for that reason that in our study we decided to combine the best of two worlds. We propose a novel hybrid optimization

41

scheme that consists of the outstanding probing capability of SA (i.e., intelligent search) with the first-rate convergence properties of TR. In our approach we first retain the results obtained over the course of multiple SA runs using the maximum number of iterations as a stopping criterion for each run (to be precise, until 50 distinct outcomes are obtained). Also to note is that other SA search parameter settings include initial temperature at the start of the algorithm in the order of 500, a fast scheduled annealing function (i.e., temperature function that is equal to initial temperature divided by k; where k denote the annealing parameter) and a re-annealing interval (number of points accepted before reannealing) that is equal to 250 [37–41]. We then employ TR using SA's results as potential starting points to carry out a search for multiple solutions. In other words, each (for a total of 50) output produced by SA is used as an initial starting point e0 and is further examined with the TR agent (instead of solely relying on the convergence properties [cooling regime] of SA [37–41]). We believe that our hybrid optimization scheme is computationally efficient compared to the ones described in [87–89] because those methods call for a higher number of starting points (obtained via a scatter-search mechanism or uniformly distributed starting points within bounds) and an exhaustive analysis using similar non-linear programming. The efficiency and effectiveness of our optimization scheme will be further discussed in Sections 4.3 and 4.4. The purpose of this hybrid optimization scheme is to solely provide insight about the relevance of the selected parameters (i.e, to provide insight regarding undesired redundancy within the model due to some irrelevant parameters [e.g., linearly dependent parameters and such]) and not to directly solve for the optimal perturbation vector. The parameters that converge to different numerical values but without significant change in the output of the objective function are highly likely to be redundant and hence some of those parameters need to be discarded in order to render the model non-redundant. In other words, a modest investigation into the flagged parameters can effectively guide the user in simplifying the model (i.e., correctly finding the DH(-) model). The suspicious parameters can be identified by simply examining the variance across the multiple solutions of a given error parameter (variance among the multiple solutions found using our SA based TR search for a given error parameter). Fig. 6 depicts a block diagram of the discussed product usage protocol. As shown in that figure, the first stage is to find the DH(-) model of the manipulator in question. The user may opt to use the analytical approach described in Section 3.1 or the automated method described in this subsection. Note that the automated approach does increase the computational burden on the user (running the SA based TR optimization scheme using available data) but is undoubtedly helpful in deriving the DH(-) model of the robot in question (in particular, it can be of great help if the geometry of the structure is complex and the guidelines of the analytical approach can no longer be extended in an obvious fashion). Hence, considering the burden of employing the automated approach in Stage 1, the analytical approach may be preferable to construct the DH(-) model, unless the robot is

Fig. 6. Top level block diagram illustrating the product usage protocol.

42

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

Table 8 Completion of Stage 1 using the automated approach: auto-relevance assessment of the various error parameters using 3DCAL and the HP20D manipulator. Selected parameters

Nominal values

Variance

a1 in mm a2 in mm a3 in mm

150.000 760.000 140.000 795.000

0.001 0.000 0.006 0.001

d4 in mm

d6 in mm δθ1 in degrees δθ 2 in degrees δθ 3 in degrees δθ 4 in degrees δθ5 in degrees δθ 6 in degrees Xtool in mm Ytool in mm Ztool in mm

105.000

1.184

0.000 0.000 0.000 0.000 0.000 0.000 2.549

2.263 0.000 0.000 0.000 0.000 2.117 0.001

0.342

0.002

70.079

1.184

kinematically complex. Provided that the DH(-) model is correctly identified using either one of those two methods, the user can complete the second stage of the calibration process; solving for the non-redundant error parameters using TR. Let us look at an example to verify the effectiveness of the automated redundance detection algorithm that we presented. We denote Table 8 as Stage 1 because it can be viewed as an intermediate procedural calibration step that the user may choose to perform at the expense of the hybrid optimization session using available data in order to deduce which of the error parameters associated to non-zero valued DH parameters are redundant. Let us also try to possibly correlate the results to the ones found using the analytical approach as the guidelines established analytically are meant to avoid redundancy as well. In other words, the goal of the two approaches (analytical and automated) is to construct the DH(-) model of a robot and hence one would expect a strong agreement among the two analyses. Table 8 presents results if one were to elect to estimate all non-zero DH parameters for the HP20D manipulator using our 3DCAL method (a total of 14 error parameters). The results were generated by employing TR as a solver using each of the 50 starting points obtained by way of SA. Other details of the experiment shall be discussed in Section 4. The objective here is to demonstrate that the established guidelines to construct the DH(-) model are in accordance with the automated assessment and hence are justifiable and well-founded. The automated evaluation (variance analysis) presented in Table 8 is consistent with the guidelines that we analytically established in Section 3.1 in order to make the model non-redundant. It follows as expected, the error corresponding to the first joint angle (δθ1) cannot be estimated using the proposed method since all calibration points are shifted by the same amount while the Euclidean distances among the points remain unaltered. This is the reason why a great level of parametric variance is observed without a change in the output of the objective function. Mathematically speaking,

∂3 y3DCAL (e) ≈ 0. ∂δθ6 ∂δd6 ∂δZtool holds true for numerous solutions of those error parameters. Based on this insight, those error parameters are kept frozen in the case of the HP20D and the relevant error parameters of the simplistic and non-redundant kinematic calibration model are identified in Stage 2 via a single run using TR (i.e., by making use of an arbitrary starting point).

4. Experimental results and discussion In this section, we present experimental results to demonstrate the effectiveness of the proposed calibration systems. We present performance analyses of 3DCAL and 1DCAL using the performance metric described in Section 2.2. The performance of our approaches are also compared to that of the state-of-the-art commercial system (MotoCal) using the same materials and performance metric. Finally, we present a discussion and thoughts for future improvements and work. 4.1. 3DCAL performance analysis In order to test our 3DCAL system, we made use of the HP20D manipulator. Let us start off by briefly describing the method we used to gather the calibration points. We used a uniformly sampled grid that is composed of 27 calibration points and that covers most of the working envelope in front of the robot. We then incorporated orientations to those 27 calibration points. Namely, orientation about the x-axis of the base frame in the order of 25 and 25°. And, orientation about the y-axis of the base frame in the order of 25 and 25°. This gave rise to a total of 135 calibration points. Fig. 7 depicts the calibration points and their orientations in the form of vectors. Unfortunately, not all calibration points could be measured using CompuGauge due to hardware limitations (occurrence of cable twist for some position measurements). This is the reason why we generated calibration points with mild orientation in the first place. After carefully testing each point, we concluded that only 87 calibration points out of the 135 are valid (because they could be measured). Table 9 shows the estimated errors using the simplified kinematic model DH(-) for the HP20D manipulator. We employed TR using an arbitrary starting point to optimize over the objective function associated with the DH(-) error model. The table shows the results for what we referred to as Stage 2 of the proposed calibration process (i.e., after removing ineffectual redundant error

∂ y3DCAL (e) ≈ 0. ∂δθ1 Also to note is the large variance in the solutions for the DH parameter d6, the last joint angle (θ6), and the Z dimension of the tool frame (Ztool). As stated in Section 3.1, the user is required to choose whether to estimate the error parameters of the last row of the DH table (e6 ) or the ones associated with tool definition (e tool ) but not both. Note that parameters d6 and the Ztool are actually linearly dependent as they describe the same dimension of the robot and tool for any configuration. That is to say

Fig. 7. Calibration points used to calibrate the HP20D manipulator [90].

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

Table 9 3DCAL and MotoCal Stage 2 calibration results using 3D CompuGauge and the DH (-) error model of the HP20D manipulator. Selected parameters

a1 in mm a2 in mm a3 in mm d4 in mm δθ 2 in degrees δθ 3 in degrees δθ 4 in degrees δθ5 in degrees

Xtool in mm Ytool in mm Ztool in mm Performance analysis Alignment test ΔX in mm ΔY in mm ΔZ in mm

Nominal values

3DCAL

MotoCal

Estimated values

Estimated values

150.000 760.000 140.000 795.000

149.901 760.971 140.071 794.966

N/A 761.000 N/A 794.960

0.000 0.000 0.000 0.000 2.549

0.135 1.001 0.156 0.614 0.046

0.143 1.002 0.156 0.614 0.066

0.342

0.015

0.020

70.079

70.615

70.625

7.4553 0.8770 1.7363

0.2798 0.3039 0.2001

0.2867 0.3117 0.2167

parameters). As mentioned in Section 3.2, we opt to estimate the tool definition in this experiment due to the fact that we did not have a well defined tool at hand. In order to obtain an initial definition (i.e., the nominal values of the tool shown in Table 9), we used the “tool calibration utility” of the DX100 controller [35]. It provides an automatic method to calculate the Xtool, Ytool and Ztool (i.e., spatial description of the mounted tool with respect to the tool flange). Obviously, it is erroneous since the robot itself was not well calibrated; and hence the need to estimate tool data. Meaning, we assume the initial estimate of the tool dimensions as nominal values and we associate the error parameters δXtool , δYtool and δZtool as described in Eq. (3) such as to refine that initial guess. The performance metric that we used to validate our approach is the so-called “alignment test” [75,79]. In order to perform this alignment test, we used eight points (the vertices of a cube that covers the majority of the dexterous workspace of the robot). Note that the chosen eight points are deemed more than sufficient to accurately evaluate the performance of a calibration scheme using the described internationally recognized standard since only 4 points are actually required provided that one of those four points is a non-coplanar testing point [75,79]. Also to note is that we restricted our approach to the dexterous workspace (as opposed to the reachable workspace [16]) because our aim is to enhance the positioning accuracy of the points that the manipulator can reach with an arbitrary orientation of the end-effector (i.e., fulfills capability required by most robotic applications with a satisfactory degree of accuracy). Those eight points are used to find the “best” fit between the coordinates of those eight points with respect to the robot's frame and measurement system's frame. The fit of the alignment procedure reflects the accuracy of the robot itself; a performance standard that is widely used in the robotics community [75,79]. Before calibration, the misfit was in the order of 7.4553, 0.770 and 1.7363 mm along the X, Y and Z dimensions, respectively. After employing our 3DCAL approach, we were able to reduce the misfit to 0.2789, 0.3039 and 0.2001 mm, respectively. According to our experimental results (shown in Table 9), our performance is comparable (in fact slightly better) compared to that of MotoCal. MotoCal has similar features. It allows user interaction (which parameters to estimate). Similar to our automated parametric relevance assessment that the user may choose to perform to

43

complete Stage 1, it also provides feedback to the user regarding which parameters ought to be kept frozen in order to have reliable results (the so-called “Acceptance Angular Stdv”) [79]. If the “Acceptance Angular Stdv” of a given parameter is above a given threshold (default value in the order of 0.04), the algorithm recommends that the parameter be omitted during the second stage of the calibration procedure. Note that this is achieved at the expense of an optimization session using available data. Based on that, in addition to θ1, DH parameters a1 and a3 are kept “frozen” during MotoCal's second calibration round. Hence, the MotoCal's results shown in Table 9 can also be referred to as Stage 2 calibration results. Note that the same data (87 calibration points) and manipulator (HP20D) were put to use to generate the results shown in Table 9 using MotoCal. It is interesting to compare the estimated values. One could state that the numerical solutions obtained using MotoCal are very similar to the ones found using our 3DCAL approach. We attribute the slight boost in performance of 3DCAL to the two extra error degrees of freedom (a1 and a3). We also present results for the other experiment where we elect to estimate the error parameters associated with the joint variables and tool data only (where we make use of yet a further simplified model; the DH(-)(-) model). Table 10 shows these results. The results are very informative for two reasons: First, it appears that most of the inaccuracy is due to errors associated with the offset of joint variables (i.e., home position). We say this because only a slight deterioration in performance is observed when comparing Tables 9 and 10. Second, it allows us to make a definite comparison between the two systems (i.e., the same number of error parameters). Based on that context, we solidly affirm that our performance is indeed comparable to that of MotoCal. To summarize the boost in performance attained using our framework with respect to MotoCal, we present Table 11. The improvements in terms of percentage are calculated using

⎛ 3DCAL Alignment Test ⎞ Improvement percent = ⎜1 − ⎟ × 100. MotoCal Alignment Test ⎠ ⎝

(15)

As stated in the latter and as it can be seen in the newly presented table, our performance is slightly better compared to that of MotoCal for the most part. The only negligible shortcoming of our 3DCAL approach is the ΔX performance result when making use of the DH(-)(-) error model of the HP20D manipulator. Table 10 3DCAL and MotoCal Stage 2 calibration results using 3D CompuGauge and the DH (-)(-) error model (home position and tool data) of the HP20D manipulator. Selected parameters

δθ 2 in degrees δθ 3 in degrees δθ 4 in degrees δθ5 in degrees Xtool in mm Ytool in mm Ztool in mm Performance analysis Alignment test ΔX in mm ΔY in mm ΔZ in mm

Nominal values

0.000 0.000 0.000 0.000 2.549

3DCAL

MotoCal

Estimated values

Estimated values

0.162 1.029 0.149 0.612 0.032

0.153 1.001 0.153 0.612 0.041

0.342

0.061

0.020

70.079

70.425

70.588

7.4553 0.8770 1.7363

0.3012 0.4624 0.3889

0.2997 0.4722 0.4039

44

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

Table 11 Summary of the improvements in performance using 3DCAL with respect to MotoCal. Alignment test

Improvement (%)

HP20D case using DH(-) model and 3D CompuGauge 2.41 ΔX 2.50 ΔY 7.66 ΔZ

In this subsection, we present similar results as in Section 4.1 but for the 1DCAL approach. The radial measurements of 105 valid calibration points were acquired using the 1D CompuGauge measurement system. We used the same procedure as for our 3DCAL system to generate the configurations (calibration points) but using the MH6 manipulator. Table 12 summarizes the estimated values and corresponding performance upon completion of Stage 2 of the calibration procedure. Note that we used a tool of known dimension in this case and hence elected to estimate the errors associated with the last link of the MH6 (e6 ) as opposed to calibrating tool data as done with the HP20D case in our 3DCAL experiment. The performance metric used here is based on the fact that given a well calibrated robot the Cartesian coordinates of a position should remain unchanged (or should change ever so slightly) as one varies the orientation of the tool frame. Hence, “Deviation” in Table 12 refers to the maximum discrepancy among the recorded radial measurements as we randomly varied the orientation of the tool for a given fixed position in space. To be precise, 34 random orientations (about all three principal axes) that are bounded by the hardware limitations of the measurement device and the robot is chosen. According to the performance standard discussed in Section 2.2, if the deviation is below 2.0 mm, it is safe to assume that the robot is fairly well calibrated for most applications [75,79]. Before calibration, a deviation in the order of 9.761 mm is reported. Based on our analytical guidelines, the offset corresponding to the first joint angle (δθ1) is the only parameter that cannot be calibrated in the case of the MH6 experiment. Thus, in this experiment, we are in favor of estimating all non-zero parameters except for θ1. Note that we completed Stage 1 for the purpose of cross validating and we are Table 12 1DCAL and MotoCal Stage 2 calibration results using 1D CompuGauge and the DH (-) error model of the MH6 manipulator.

a1 in mm a2 in mm a3 in mm

in degrees in degrees in degrees in degrees in degrees

Performance analysis Deviation in mm

4.2. 1DCAL performance analysis

Nominal values

Selected parameters

δθ 2 δθ 3 δθ 4 δθ5 δθ 6

HP20D case using DH(-)(-) model and 3D CompuGauge 0.5 ΔX 2.07 ΔY 3.71 ΔZ

Selected parameters

Table 13 1DCAL and MotoCal Stage 2 calibration results using 1D CompuGauge and the DH (-)(-) error model (home position) of the MH6 manipulator.

1DCAL

MotoCal

Estimated values

Estimated values

148.140 614.381 155.281 640.534

N/A 614.748 N/A N/A 94.827

d4 in mm

150.000 614.000 155.000 640.000

d6 in mm δθ 2 in degrees δθ 3 in degrees δθ 4 in degrees δθ5 in degrees δθ 6 in degrees

95.000

94.855

0.000 0.000 0.000 0.000 0.000

0.128 0.046 0.117 0.915 0.837

Performance analysis Deviation in mm

9.761

1.131

0.199 0.017 0.136 0.906 0.909

1.301

Nominal values

0.000 0.000 0.000 0.000 0.000

9.761

1DCAL

MotoCal

Estimated values

Estimated values

0.188 0.083 0.065 0.914 0.882

1.372

0.353 0.074 0.125 0.907 0.855

1.434

pleased to make note that the auto-assessment does again conform with the analytical approach. After executing Stage 2 (calibration using the DH(-) model of the MH6 manipulator), we were able to reduce the deviation to 1.131 mm using our 1DCAL approach. The commercial system MotoCal can also calibrate a given robot with the supply of radial measurements [79]. Note that MotoCal's algorithm has recommended not to consider 3 of the non-zero parameters shown in Table 12. Unlike the 3D case, our numerical results greatly differ compared to the ones generated using MotoCal. This could be due to the extra 3 parameters in the case of 1DCAL. Note that in this particular scenario, our approach has outperformed MotoCal. To make a definite comparison to our best ability, we present Table 13 where only the joint angle offsets (i.e., calibration using the DH(-)(-) error model) are estimated. The results from both systems remain different but a slight increase in performance is attained when employing 1DCAL. Additionally, by comparing the performance reported in the two tables (Tables 12 and 13), we re-affirm that indeed most positioning inaccuracy is mostly due to incorrect joint variable offset values (i.e., de-calibrated home position of the robot) because only a slight decline in performance is noted. We conclude this subsection by presenting a similar examination like that of Table 11 but with our 1D approach. Table 14 clearly shows that 1DCAL yields improved results in both cases. 4.3. Discussion A well known issue that needs to be discussed when examining the performance of a calibration scheme is the sensitivity of the system to the number of calibration points. After experimenting, we found that as long as there are at least twice as many calibration points as calibration parameters (error degrees of freedom) both systems (3DCAL and 1DCAL using the recommended solver over the course of the two calibration stages) converged to the same solution. This does not hold true with MotoCal. Meaning, MotoCal requires a comparatively greater number of calibration points in order to yield the same calibration results. It is instructive to note that one needs to incorporate at least two different wrist orientations for each calibration point at different spatial positions. Also to note is that, the configurations used to calibrate the robot must not be singular (since the corresponding calibration points may be mapped to infinitely many solutions or the manipulator's ability to adequately and effectively Table 14 Summary of the improvements in performance using 1DCAL with respect to MotoCal. Performance metric and error models

Improvement (%)

Deviation using DH(-) model Deviation using DH(-)(-) model

13.07 4.32

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

position and orient a tool has been compromised [12,15–19]). As promised to the reader, here we shall discuss how to efficiently define the solution space Ω for the 3DCAL and 1DCAL approaches. A manipulator requires calibration for two main reasons: (1) if the accuracy of the robot (out of the box) is deemed inadequate (further improvement in accuracy is required by the application) or (2) if the robot has been repaired (after collision or parts replacement is performed). Aforementioned, the dimensions of the solution space Ω are manufacturer specific. For example, Yaskawa, Motoman Robotics Inc., manipulators are known to be fairly accurate “out of the box.” Hence, if further improvement in accuracy is required, a conservative restriction on the solution space is sufficient. To be specific, the bounds related to angular error parameters, − e and e¯ (such as home position), may be limited to 71.0° and the ones corresponding to length error parameters (such as link lengths and link offsets) may be limited to 71.0 mm. Robot manufacturers typically also provide a service for re-building and repairing their manipulators when needed (the second reason discussed in the latter). In this second scenario, it is first recommended that the user make use of the “pin alignment” technique to visually (with the “naked eye”) get a rough estimate of the home position. Then, the user may employ the presented calibration algorithms subject to a comparatively wider solution space. That is, for this particular scenario, we recommend that the bounds that corresponds to angular error parameters, − e and e¯ , ought to be limited to 73.0° and the ones bounding link length/offset error parameters may be limited to 72.0 mm. It is important to note that our recommendation here on how to set Ω are based on Yaskawa, Motoman Robotics, Inc. manufacturing and repairing capabilities. Thus, those bounds need to be revised accordingly as they are manipulator manufacturer and servicing quality dependent. In addition to the specification of the solution space Ω , our 1DCAL approach requires the definition of the spherical volumetric solution space V. Note that the discrepancy between the set of calculated values r and the set of measured values r˜ is due to decalibration. In fact, both, Ω and V, are functions of the level of decalibration. That is the reason why we recommend to set Ω in accordance with the two scenarios discussed in the above. We recommend a more methodological approach to estimate the radius of V. That is, the spherical volumetric solution space V can be readily defined by computing its radius as

rV = 2 × max (r − r˜ ), rV ∈

where r and r˜ are sets of calculated and measured values. We also propose that V be centered at p⁎e0 . V must include the location of the encoder and hence we believe it is wise to incorporate a safety factor (in particular since p⁎e0 is an initial guess). The largest measurement discrepancy is multiplied by 2 in order to ensure that the actual position of the encoder is within the volumetric space V. To compare the degree to which our hybrid optimization scheme (employed in Stage 1 if the user chooses to make use of the automated approach) is successful in providing useful relevance assessment of the selected error parameters effectively, we took on the task of implementing a similar optimization scheme. Namely, we looked at the “Multi-start” framework described in [87]. We made use of the scatter search method to generate potential starting points and employed TR for convergence purposes (i.e., we used different methods [SA versus Multi-start] to generate the starting points but employed an identical solver to analyze each initial condition). In order for that framework (Multi-start based TR) to convey a comparable evaluation (i.e., as indicative), a scatter search composed of over 450 starting points is required. For example, referring back to Table 8, our technique (SA based TR approach using 50 initial conditions) indicated a variance in the order of 2.263 for the error parameter

45

δθ1 upon completion of the first stage of calibration. Generating the same amount of starting points using the method described in [87] and performing a relevance evaluation on that parameter returns a variance analysis in the order of 0.981 only. Note that it also took a slightly higher number of iterations to converge for a given starting point (additional 7–9 iterations in some cases). Our framework is obviously weakly dependent on that particular parameter for the reasons explained in Section 3.1 (i.e., a change in that parameter has null impact in the output of the objective function). We attribute our higher ranking performance in this domain (i.e., relevance evaluation and convergence rate) to SA's ability not only to explore the solution space V thoroughly but also to output starting points that are closer to the minima of the cost function. To characterize the objective function associated with the DH(-) error model, we performed the automated saliency assessment on its error parameters. We made note of the following interesting observations. All 50 potential starting points converged to the same numerical solution for the bounded search space Ω (i.e., the automated evaluation indicates no variance). Also to note is that identical results were seen when employing the Multi-start framework described in [87] (instead of SA) to generate a different set of “casual” starting points that is significantly higher in number. This empirical study proves that after the completion of the first stage, almost all degeneracy in the error model has been eliminated and the cost function associated with that model is well-conditioned. In other words, the objective function exhibits a convex like characteristic (within Ω ) [91]. Hence, the DH(-) error models warrant robust calibration results as this check verifies that the obtained multiple numerical solutions are consistent regardless the initial position of a widespread starting points. This is the reason why we recommend a single run using TR with an arbitrary starting point in the second stage of the accuracy improvement procedure as convergence to a reliable and robust solution is inevitable (i.e., a SA based TR search for multiple minima is simply unnecessary). 4.4. Future work A valuable area of future work would be to improve processing speed. The current prototype MATLAB implementation, using no hardware acceleration of any kind, consumes approximately on average 271 section if the user opts to execute Stage 1 of the calibration procedure (employing SA based TR optimization scheme). Note that making use of the “Multi-start” approach discussed in Section 4.3 to complete this stage and attain a comparable performance in characterizing the error parameters consumes approximately on average 3018 section (since that approach requires employing TR using 450 “casual” starting points to yield a comparable performance in that domain). Similar to the improvement in accuracy that we calculated using Eq. (15), the improvement in efficiency (in terms of time consumption) using our novel optimization scheme with respect to the “Multi-Start” based TR framework can be calculated. Our SA based TR approach is found to be 91.02% more efficient with respect to the “Multi-Start” based TR framework. Stage 2 of the calibration procedure using DH(-) and DH(-)(-) models consumes on average 7.6 and 6.8 s, respectively, using TR and a single arbitrary starting point. The work is implemented on an HP Personal Computer (PC) with processor speed of 3.22 GHz. The obvious thing to note is that the reported temporal rates are subject to the amount of calibration-points used to calibrate the robot in question. Also to note is that the calculated time consumption appears to be less on average in the case of the 1DCAL calibration system (in comparison with 3DCAL). In order to accelerate the automated approach of Stage 1, we would like to refer the reader to [92]. In that work, we have proven that as long as processes are parallel, they can be greatly accelerated. The process related to each starting point is independent from another

46

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

(i.e., parallel); and hence, those processes can be simultaneously done on a platform like the one described in [92]. We are also in the process of exploring higher quality metrology equipments (compared to CompuGauge) in order to evaluate the peak performance and limitations of the presented frameworks. We also believe that extending the presented simplistic approach so as to include other sources of error is another potentially fruitful area of future research (i.e., to devise methodologies to construct non-redundant error models that take into account parameters related to elasticity of joints and links, load proportion and joint deflection, etc.). It is also our hope to investigate computationally efficient kinematic and dynamic related calibration schemes for redundant manipulators (i.e., robots that possess more than six degrees of freedom) like the Yaskawa Motoman Robotics, Inc. IA20 that possess 7 revolute joints [93].

parameters and construct the DH(-) error model for a given manipulator. We have conducted experiments that show that the automated evaluation (variance analysis) is consistent with the guidelines that we analytically established in order to identify irrelevant parameters and to render the model non-redundant.

Acknowledgments The authors would like to thank Mr. Wade Hickle of Yaskawa Motoman Robotics, Inc. for providing the materials used in this study. We also wish to express our appreciation to personnel of Yaskawa Motoman Robotics, Inc. for helping with the experiments. Finally, we wish to thank the editor and the anonymous reviewers for helping to strengthen this paper.

5. Conclusion In this paper, we present new computational efficient kinematic calibration algorithms for industrial robots. These include the 3DCAL system that uses Cartesian coordinates and the 1DCAL approach that requires the supply of radial measurements only. Both methods make use of an error model that is based on the original DH convention. Both systems make use of a two stage procedure. The first stage consists of properly deducing the parameterically continuous and non-redundant DH(-) kinematic calibration model of a manipulator. That can be achieved either using an analytical or an automated approach. Note that the automated approach does increase the computational burden on the user but is very helpful in deriving the DH(-) model of the robot in question. In particular, it can greatly help if the geometry of the structure is complex and cannot be easily interpreted analytically. After completing Stage 1 (i.e., constructing a redundancy free and continuous error model), TR is brought into play to carry out the second stage of the calibration process. Based on our experiment, we conclude that most of positioning inaccuracy of manipulators like the HP20D and MH6 is due to the error parameters related to joint angles (i.e., home position of robot). Hence, we have also introduced the DH(-)(-) error model where the home position (and tool if needed) only is calibrated. Note that those two models are less complex (compared to the numerous error models discussed in Section 1) and hence do not require burdensome computation. Although incomplete, in this work, we have demonstrated that those simplified error models (DH(-) and DH(-)(-)) in conjunction with the proposed optimization scheme (TR) are capable of exceptional performance comparable or superior to that of an existing commercial system. Using our 3DCAL, we are able to successfully calibrate the Yaskawa Motoman Robotics, Inc. 6DOF serial manipulator HP20D to a satisfactory degree for most applications (according to performance standard and metric described in Section 2.2). Similarly, using the presented 1DCAL approach, we have also adequately calibrated the MH6 robot. We make use of a novel hybrid optimization scheme in our automated approach (Stage 1 of calibration procedure) to effectively target conspicuous error parameters. The scheme consists of first collating a rich set of potential starting points using SA. We then employ the solver TR using those starting points to carry out a search for multiple optimizers. To our best knowledge, we have not seen that done in any other work. We found our optimization technique not only to be beneficial from a computational stand point of view (compared to other search mechanism that call for a higher number of starting points) but also to be very well suited for the first stage in this application as it can provide an automated parametric relevance assessment to the end-user. Based on that insight, the user can easily identify the ineffectual/redundant error

Appendix A. Expressional details of homogeneous transformation matrices

⎡1 0 0 0⎤ ⎥ ⎢ 0 cos sin 0⎥ ( ) − ( ) α α i i Rot (x, αi ) = ⎢ Rot (z, θ i ) ⎢0 sin (αi ) cos (αi ) 0⎥ ⎥ ⎢ ⎣0 0 0 1⎦ ⎡cos (θ i ) − sin (θ i ) 0 0⎤ ⎥ ⎢ sin (θ i ) cos (θ i ) 0 0⎥ Trans (z, di ) =⎢ ⎢ 0 0 1 0⎥⎥ ⎢ ⎣ 0 0 0 1⎦ ⎡1 0 0 0 ⎤ ⎥ ⎢ 0 1 0 0⎥ Trans (x, ai ) =⎢ ⎢ 0 0 1 di ⎥ ⎥ ⎢ ⎣0 0 0 1 ⎦ ⎡1 ⎢ = ⎢0 ⎢0 ⎢ ⎣0

0 0 ai ⎤ ⎥ 1 0 0 ⎥ A (q , e ) i i i 0 1 0⎥ ⎥ 0 0 1⎦

= Rot (z, θ i + δθ i ) Trans (z, di + δdi ) Trans (x, ai + δai ) Rot (x, αi + δαi ) A i ( qi , e i ) = ⎡ c (θ i + δθ i ) − s (θ i + δθ i ) c (α i + δα i ) s (θ i + δθ i ) s (α i + δα i ) (a i + δa i ) c (θ i + δθ i ) ⎤ ⎢ ⎥ ⎢ s (θ i + δθ i ) c (θ i + δθ i ) c (α i + δα i ) − c (θ i + δθ i ) s (α i + δα i ) (a i + δa i ) s (θ i + δθ i ) ⎥ ⎢ ⎥ 0 s (α i + δα i ) c (α i + δα i ) d i + δd i ⎢ ⎥ ⎢⎣ ⎥⎦ 0 0 0 1

RPY (RXtool, RYtool, RZtool ) ⎡ c (RZtool ) s (RXtool ) s (RYtool ) ⎢c (RZtool ) c (RYtool ) ⎢ − c (RXtool ) s (RZtool ) ⎢ ⎢ c (RYtool ) s (RZtool ) c (RZtool ) c (RXtool ) + s (RZtool ) ⎢ s (RXtool ) s (RYtool ) ⎢ ⎢ c (RYtool ) s (RXtool ) ⎢ − s (RYtool ) ⎢⎣ 0 0

⋯

⋯

⎤ ⎥ s (RZtool ) s (RXtool ) + c (RZtool ) c (RXtool ) s (RYtool ) 0⎥ ⎥ c (RXtool ) s (RZtool ) s (RYtool ) − c (RZtool ) s (RXtool ) 0⎥ ⎥ c (RXtool ) c (RYtool ) 0⎥ 0 1⎥ ⎥ ⎥⎦ where c () and s () are the short hands of cos() and sin(), respectively.

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

References [1] S.Y. Nof, Handbook of Industrial Robotics, vol. 1, John Wiley & Sons, Inc., New York, NY, 1999. [2] B.W. Mooring, Z.S. Roth, M.R. Driels, Fundamentals of Manipulator Calibration, Wiley, New York, 1991. [3] A. Nubiola, I.A. Bonev, Absolute calibration of an abb irb 1600 robot using a laser tracker, Robot. Comput.-Integr. Manuf. 29 (2013) 236–245. [4] B. Horn, Robot Vision, MIT Press, Cambridge, Massachusetts and London, England, 1986. [5] P.G. Mikell, R. Mitchell Weis, N.G.O.N. Nagel, Industrial Robotics Technology, Programming and Applications, 1986. [6] A. Nubiola, I.A. Bonev, Absolute robot calibration with a single telescoping ballbar, Precis. Eng. (2014). [7] B. Karan, M. Vukobratović, Calibration and accuracy of manipulation robot models an overview, Mech. Mach. Theory 29 (1994) 479–500. [8] Y. Andrew Liou, P.P. Lin, R.R. Lindeke, H.-D. Chiang, Tolerance specification of robot kinematic parameters using an experimental design technique the Taguchi method, Robot. Comput.-Integr. Manuf. 10 (1993) 199–207. [9] C. Gong, J. Yuan, J. Ni, Nongeometric error identification and compensation for robotic system by inverse calibration, Int. J. Mach. Tools Manuf. 40 (2000) 2119–2137. [10] M. Abderrahim, A. Khamis, S. Garrido, L. Moreno, Accuracy and calibration issues of industrial manipulators, Ind. Robot.: Program. Simul. Appl. (2006). [11] A. Elatta, L.P. Gen, F.L. Zhi, Y. Daoyuan, L. Fei, An overview of robot calibration, Inf. Technol. J. 3 (2004) 74–78. [12] L. Sciavicco, B. Siciliano, Modelling and Control of Robot Manipulators, Springer-Verlag, London, 2000. [13] M.R. Driels, U.S. Pathre, Robot calibration using an automatic theodolite, Int. J. Adv. Manuf. Technol. 9 (1994) 114–125. [14] J.M.S. Motta, G.C. de Carvalho, R. McMaster, Robot calibration using a 3d vision-based measurement system with a single camera, Robot. Comput.-Integr. Manuf. 17 (2001) 487–497. [15] R.P. Paul, Robot Manipulators: Mathematics, Programming, and Control: The Computer Control of Robot Manipulators, 1981. [16] M.W. Spong, S. Hutchinson, M. Vidyasagar, Robot Modeling and Control, John Wiley & Sons, New York, 2006. [17] J. Robert, Fundamentals of Robotics: Analysis and Control, Prentice-Hall, Canada, 1990. [18] J.J. Craig, Introduction to Robotics, vol. 7, Addison-Wesley Reading, MA, 1989. [19] M. Brady, Robot Motion: Planning and Control, MIT Press, Cambridge, Massachusetts and London, England, 1982. [20] M. Abderrahim, A. Whittaker, Kinematic model identification of industrial manipulators, Robot. Comput.-Integr. Manuf. 16 (2000) 1–8. [21] X. Yang, L. Wu, J. Li, K. Chen, A minimal kinematic model for serial robot calibration using poe formula, Robot. Comput.-Integr. Manuf. 30 (2014) 326–334. [22] J. Denavit, A kinematic notation for lower-pair mechanisms based on matrices, Trans. ASME. J. Appl. Mech. 22 (1955) 215–221. [23] D. Whitney, C. Lozinski, J.M. Rourke, Industrial robot forward calibration method and results, J. Dyn. Syst. Meas. Control 108 (1986) 1–8. [24] H.W. Stone, A.C. Sanderson, A prototype arm signature identification system, in: 1987 IEEE International Conference on Robotics and Automation. Proceedings, vol. 4, IEEE, Raleigh, N.C, 1987, pp. 175–182. [25] S.A. Hayati, Robot arm geometric link parameter estimation, in: The 22nd IEEE Conference on Decision and Control, 1983, vol. 22, IEEE, San Antonio, TX, 1983, pp. 1477–1483. [26] S. Hayati, M. Mirmirani, Improving the absolute positioning accuracy of robot manipulators, J. Robot. Syst. 2 (1985) 397–413. [27] G. Du, P. Zhang, Online robot calibration based on vision measurement, Robot. Comput.-Integr. Manuf. 29 (2013) 484–492. [28] H.W. Stone, Kinematic Modeling, Identification, and Control of Robotic Manipulators, Kluwer Academic, Boston, 1987. [29] W. Veitschegger, C.-H. Wu, Robot accuracy analysis based on kinematics, IEEE J. Robot. Autom. 2 (1986) 171–179. [30] H. Zhuang, Z.S. Roth, F. Hamano, A complete and parametrically continuous kinematic model for robot manipulators, IEEE Trans. Robot. Autom. 8 (1992) 451–463. [31] K. Okamura, F. Park, Kinematic calibration using the product of exponentials formula, Robotica 14 (1996) 415–422. [32] R. He, Y. Zhao, S. Yang, S. Yang, Kinematic-parameter identification for serialrobot calibration based on poe formula, IEEE Trans. Robot. 26 (2010) 411–423. [33] X. Yang, L. Wu, J. Li, K. Chen, A minimal kinematic model for serial robot calibration using poe formula, Robot. Comput.-Integr. Manuf. 30 (2014) 326–334. [34] Yaskawa Motoman Robotics, Inc., “NX100 Controller Manual”, Part Number:149201-1, Yaskawa Motoman Robotics Inc., Yaskawa America, Inc., Motoman Robotics Division, Headquarters and Manufacturing Facility, 100 Automation Way, Miamisburg, OH 45342, United States, 2007. URL: 〈www.moto man.com/〉. [35] Yaskawa Motoman Robotics, Inc., “DX100 Controller Manual”, Part Number:156431-1CD, Yaskawa Motoman Robotics Inc., Yaskawa America, Inc., Motoman Robotics Division, Headquarters and Manufacturing Facility, 100 Automation Way, Miamisburg, OH 45342, United States, 2011. URL: 〈www. motoman.com/〉.

47

[36] Yaskawa Motoman Robotics, Inc., “FS100 Controller Manual”, Part Number:159550-1CD, Yaskawa Motoman Robotics Inc., Yaskawa America, Inc., Motoman Robotics Division, Headquarters and Manufacturing Facility, 100 Automation Way, Miamisburg, OH 45342, United States, 2011. URL: 〈www. motoman.com/〉. [37] S. Kirkpatrick, Optimization by simulated annealing: quantitative studies, J. Stat. Phys. 34 (1984) 975–986. [38] L. Davis, Genetic Algorithms and Simulated Annealing, 1987. [39] P.J.M. Van Laarhoven, E.H. Aarts, Simulated Annealing: Theory and Applications, D. Reidel Publishing Company, Dordrecht, 1987. [40] E. Aarts, J. Korst, Simulated Annealing and Boltzmann Machines: A Stochastic Approach to Combinatorial Optimization and Neural Computing, 1988. [41] T.R. Tuinstra, Automatic segmentation of small pulmonary nodules in computed tomography data using a radial basis function neural network with application to volume estimation (Ph.D. thesis), University of Dayton, 2008. [42] J.J. Moré, D.C. Sorensen, Computing a trust region step, SIAM J. Sci. Stat. Comput. 4 (1983) 553–572. [43] R.H. Byrd, R.B. Schnabel, G.A. Shultz, Approximate solution of the trust region problem by minimization over two-dimensional subspaces, Math. Program. 40 (1988) 247–263. [44] T.F. Coleman, Y. Li, An interior trust region approach for nonlinear minimization subject to bounds, SIAM J. Optim. 6 (1996) 418–445. [45] T. Steihaug, The conjugate gradient method and trust regions in large scale optimization, SIAM J. Numer. Anal. 20 (1983) 626–637. [46] K.M. Passino, Biomimicry for Optimization, Control, and Automation, Springer-Verlag, London, 2005. [47] Yaskawa Motoman Robotics, Inc., HP20D Series: Technical Specifications, Yaskawa Motoman Robotics Inc., 2012. URL: 〈www.motoman.com/datasheets/ HP20D_HP20D-6.pdf/〉. [48] Yaskawa Motoman Robotics, Inc., MH6 Series: Technical Specifications, Yaskawa Motoman Robotics Inc., 2011. URL: 〈www.motoman.com/datasheets/ MH6_MH6S.pdf/〉. [49] Yaskawa Motoman Robotics, Inc., Motoman MH5: General Purpose and Handling, Yaskawa Motoman Robotics Inc., 2011. URL: 〈www.motoman.co.uk/ uploads/tx_catalogrobot/MH5_E.pdf/ 〉. [50] K. Lau, R. Hocken, W. Haight, Automatic laser tracking interferometer system for robot metrology, Prec. Eng. 8 (1986) 3–8. [51] A. Nubiola, M. Slamani, I.A. Bonev, A new method for measuring a large set of poses with a single telescoping ballbar, Prec. Eng. 37 (2013) 451–460. [52] B. Mooring, S. Padavala, The effect of kinematic model complexity on manipulator accuracy, in: 1989 IEEE International Conference on Robotics and Automation, 1989. Proceedings, IEEE, , Arizona, USA. 1989, pp. 593–598. [53] G. Puskorius, L. Feldkamp, Global calibration of a robot/vision system, in: 1987 IEEE International Conference on Robotics and Automation. Proceedings, vol. 4, IEEE, Raleigh, N.C., 1987, pp. 190–195. [54] Y. Meng, H. Zhuang, Self-calibration of camera-equipped robot manipulators, Int. J. Robot. Res. 20 (2001) 909–921. [55] C.S. Gatla, R. Lumia, J. Wood, G. Starr, An automated method to calibrate industrial robots using a virtual closed kinematic chain, IEEE Trans. Robot. 23 (2007) 1105–1116. [56] F. Boochs, R. Schutze, C. Simon, F. Marzani, H. Wirth, J. Meier, Increasing the accuracy of untaught robot positions by means of a multi-camera system, in: 2010 International Conference on Indoor Positioning and Indoor Navigation (IPIN), IEEE, Zurich, Switzerland, 2010, pp. 1–9. [57] J.H. Jang, S.H. Kim, Y.K. Kwak, Calibration of geometric and non-geometric errors of an industrial robot, Robotica 19 (2001) 311–321. [58] R.P. Judd, A.B. Knasinski, A technique to calibrate industrial robots with experimental verification, IEEE Trans. Robot. Autom. 6 (1990) 20–30. [59] J. Chen, L.-M. Chao, Positioning error analysis for robot manipulators with all rotary joints, IEEE J. Robot. Autom. 3 (1987) 539–545. [60] J.F. Jarvis, Microsurveying: towards robot accuracy, in: 1987 IEEE International Conference on Robotics and Automation. Proceedings, vol. 4, IEEE, Raleigh, N. C., 1987, pp. 1660–1665. [61] H.W. Stone, A.C. Sanderson, C.P. Neuman, Arm signature identification, in: 1986 IEEE International Conference on Robotics and Automation. Proceedings, vol. 3, IEEE, San Francisco, CA, 1986, pp. 41–48. [62] W. Khalil, S. Besnard, Geometric calibration of robots with flexible joints and links, J. Intell. Robot. Syst. 34 (2002) 357–379. [63] G. Schiele, R. Hofmann, T.v. Diep, H. van Kunzmann, F. Wäldele, K. Busch, Transportable measuring instrument for testing the positional accuracy of a program-controlled appliance arm, German Patent No. DE 3504464, 1986, C1. [64] N. Vira, K. Lau, An extensible ball bar for evaluation of robots' positioning performance, J. Robot. Syst. 4 (1987) 799–814. [65] M. Driels, Using passive end-point motion constraints to calibrate robot manipulators, J. Dyn. Syst. Meas. Control 115 (1993) 560–566. [66] N. Juneja, A.A. Goldenberg, Kinematic calibration of a re-configurable robot (robotwin), in: 1997 IEEE International Conference on Robotics and Automation, 1997. Proceedings, vol. 4, IEEE, Albuquerque, New Mexico, 1997, pp. 3178– 3183. [67] H. Ota, T. Shibukawa, T. Tooyama, M. Uchiyama, Forward kinematic calibration and gravity compensation for parallel-mechanism-based machine tools, Proc. Inst. Mech. Eng. Part K: J. Multi-body Dyn. 216 (2002) 39–49. [68] L. Beyer, J. Wulfsberg, Practical robot calibration with ROSY, Robotica 22 (2004) 505–512. [69] H. Zhuang, L.K. Wang, Z.S. Roth, Error-model-based robot calibration using a modified cpc model, Robot. Comput.-Integr. Manuf. 10 (1993) 287–299.

48

T. Messay et al. / Robotics and Computer-Integrated Manufacturing 37 (2016) 33–48

[70] A. Nubiola, M. Slamani, A. Joubair, I.A. Bonev, Comparison of two calibration methods for a small industrial robot based on an optical cmm and a laser tracker, Robotica (2013) 1–20. [71] A. Goswami, A. Quaid, M. Peshkin, Complete parameter identification of a robot from partial pose information, in: 1993 IEEE International Conference on Robotics and Automation, 1993. Proceedings, IEEE, Atlanta, Georgia, 1993, pp. 168–173. [72] L. Foulloy, R.B. Kelley, Improving the precision of a robot, in: 1984 IEEE International Conference on Robotics and Automation. Proceedings, vol. 1, IEEE, Atlanta, Georgia, 1984, pp. 62–67. [73] W. Veitschegger, C.-H. Wu, A method for calibrating and compensating robot kinematic errors, in: 1987 IEEE International Conference on Robotics and Automation. Proceedings, vol. 4, IEEE, Raleigh, N.C, 1987, pp. 39–44. [74] T.-W. Hsu, L.J. Everett, Identification of the kinematic parameters of a robot manipulator for positional accuracy improvement, in: Proceedings of the Computation in Engineering Conference, 1985, pp. 263–267. [75] Dynalog Inc., 3-D COMPUGAUGE Manual (version 4.0), Compugauge: Dynalog Inc., Dynalog, Inc. 6001 1 N. Adams Road, Suite 125, Bloomfield Hills, MI 48304, United States, 2002. URL: 〈www.dynalog-us.com/ 〉. [76] K. Levenberg, A method for the solution of certain problems in least squares, Q. Appl. Math. 2 (1944) 164–168. [77] D.W. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, J. Soc. Ind. Appl. Math. 11 (1963) 431–441. [78] J.J. Moré, The Levenberg–Marquardt algorithm: implementation and theory, in: Numerical Analysis, Springer-Verlag, Berlin, 1978, pp. 105–116. [79] Yaskawa Motoman Robotics, Inc., MotoSoft:“MotoCal 3.0.1 User's Manual”, Part Number:140724-1, Revision:2, Yaskawa Motoman Robotics Inc., Yaskawa America, Inc., Motoman Robotics Division, Headquarters and Manufacturing Facility, 100 Automation Way, Miamisburg, OH 45342, United States, 2007. URL: 〈www.motoman.com/〉. [80] H. Zhuang, Self-calibration of parallel mechanisms with a case study on stewart platforms, IEEE Trans. Robot. Autom. 13 (1997) 387–397. [81] J. Gentle, Matrix Algebra: Theory, Computations, and Applications in Statistics, 2007.

[82] A. Jeffrey, Advanced Engineering Mathematics, Harcourt/Academic Press, Burlington, Massachusetts, 2001. [83] R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, New York, NY, 2012. [84] D. Lay, Linear Algebra and its Applications, Pearson Education, Inc., New York, NY, 2003. [85] F.V. Berghen, Levenberg–Marquardt algorithms vs trust region algorithms, IRIDIA, Université Libre de Bruxelles, 2004. [86] G. Alici, B. Shirinzadeh, Enhanced stiffness modeling, identification and characterization for robot manipulators, IEEE Trans. Robot. 21 (2005) 554–564. [87] Z. Ugray, L. Lasdon, J. Plummer, F. Glover, J. Kelly, R. Martí, Scatter search and local nlp solvers: a multistart framework for global optimization, INFORMS J. Comput. 19 (2007) 328–340. [88] F. Glover, A template for scatter search and path relinking, in: Artificial Evolution, Springer, Berlin, Germany, 1998, pp. 1–51. [89] L. Dixon, G. Szegö, The global optimization problem: an introduction, Towards Global Optim. 2 (1978) 1–15. [90] Y.E. Corp., S. Corp., “MotoSimEG”, Version 3.50, Yaskawa Motoman Robotics Inc., Yaskawa America, Inc., Motoman Robotics Division, Headquarters and Manufacturing Facility, 100 Automation Way, Miamisburg, OH 45342, United States, 2009. URL: 〈www.motoman.com/ 〉. [91] S.P. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, New York, NY, 2004. [92] T. Messay, C. Chen, R. Ordóñez, T.M. Taha, Gpgpu acceleration of a novel calibration method for industrial robots, in: Proceedings of the 2011 IEEE National Aerospace and Electronics Conference (NAECON), IEEE, Dayton, Ohio, 2011, pp. 124–129. [93] Yaskawa Motoman Robotics, Inc., IA20 Manipulator Manual., Yaskawa Motoman Robotics Inc., 2007. URL: 〈www.motoman.com/motomedia/manuals/ docs_151529-1.pdf/〉.