Development and Experimental Studies of an Identification Method of Kinematic Parameters for Industrial Robots without External Measuring Instruments

This paper proposes a method for the identification of kinematic parameters of multilink series industrial robots, which does not require the use of complex and expensive equipment for high-precision external measurements of the position and orientation of an end effector in a Cartesian coordinate system. This method, by means of simple and affordable tools, enables us to achieve a significantly increased accuracy of movement of end effectors in serial robots performing various technological operations. The proposed method is experimentally verified and can be applied directly in production lines.


Introduction
In the present day, the conventional approach to industrial robots (IR) programming is manual end effector (EE) driving to required points in a working area with the help of a teach pendant. These points become base points of the EE movement trajectory considering actual location of machined parts and other additional equipment. An IR controller saves vectors of generalized coordinates (rotation angles or linear displacements of IR joints), which correspond to these base points of trajectory. Then, depending on the problem being solved, with help of various types of interpolation, the required trajectories of IR movement passing through the specified base points are formed. Accuracy of EE movement along these trajectories in automatic mode is defined by the accuracy of all actuators, which track desired values of generalized coordinates. This problem is well known, and, today, there are a lot of methods which can provide high-accuracy control of IR actuators [1][2][3][4][5][6].
In modern applications, IR preforms operations in uncertain working environments and generates EE trajectories offline, or by means of information from their sensors (computer vision systems, probes and other). In this case, the base points of the trajectory are set in a base coordinate system (BCS) of IR. Real positioning accuracy of EE is defined by the accuracy of the IR kinematic model, since the IR controller uses it for calculation of EE position, by means of rotation angle information or linear displacements of IR joints. This is why, if the parameters of the IR kinematic model are used in the controller differ from their real values, then the difference between the calculated EE position and its real position in BCS can be several millimeters. This is unacceptable for many technological operations [7][8][9][10][11][12][13][14].
These methods have the following disadvantages: poor usability in real production conditions, requiring specially trained employees, requiring information about hardlyobtainable parameters, for example, a covariation matrix, and the high cost of measurement equipment. This equipment is often inaccessible for system integrators. As a result, there is the task of developing simplified procedures for the identification of IR kinematic parameters that does not require the use of expensive specialized measurement equipment.
In [36], an identification method of geometric parameters for multilink manipulators is offered. The identification process is performed by means of EE driving to the same point. This point is a small hole on a metallic sheet. During experimental research, this point can be reached by EE with different orientations up to 40 times. This allows us to form 20 pairs of points and write 60 equations for the identification of IR parameters. As a result, the improvement parameters allow us to decrease the EE positioning error from 1 sm to 1 mm. However, experience has shown that to decrease this error to 0.1-0.3 mm, it is necessary to use about 600 pairs of points. To get such s number of points by means of a specified approach is difficult enough. Probably for this reason, this method has not been popular.
In [37], data for the identification of IR kinematic parameters are obtained in the following manner. A special probe with known parameters is placed into holes in a calibrating plate. Requirements for the production of this probe are very high as it needs to know the position of all hole center points precisely for the forming of equations. Also, touch sensors can be used to define the contact of the probe with the hole bottom. There were 102 measurements obtained during experimental research. After performing the identification procedure and applying updated parameters in the IR kinematic model, the deviation of EE from the desired position was decreased to 0.3 mm. In the same manner, identification of kinematic parameters for UR series robots is performed.
In [38], a method for the identification of IR kinematic parameters using a contact probe and three spheres mounted on a platform capable of changing its orientation is proposed. During the identification process, the probe is driven into 15 points of each sphere for each of five fixed orientations of platform. Data for the angles of rotation of the hinges, information about the diameter of the spheres, and the distances between them are used to adjust IR kinematic parameters. This method does not require the use of external measuring devices, but its disadvantage is the need to use additional precisely manufactured equipment.
In [39], a method is proposed for the identification of parameters for manipulators, based on minimizing the deviation of distances between points of trajectories specified in the working area, and corresponding distances were calculated using its kinematic model. At the same time, ref. [40] describes a variant of this approach using only one point of a manipulator's working area with known coordinates. A disadvantage of this approach, one needs one or more points in the workspace, with known coordinates in the manipulator BCS. However, this is not always possible since it involves the use of external measuring devices and special procedures for setting these points. The disadvantages also include the empirical selection of correction coefficients.
In [17,41,42], similar methods are offered. In these methods, data for the identification procedure are formed by means of reaching EE to points of the same plane. In this case, data collection for identification is simpler, and can be automated. However, the result of such identification is worse than previous cases. It can be explained that the equation of the plane has four unknown parameters, and when EE reaches the point of this plane, one can only get three coordinates. To solve this problem, we must make some assumptions that may worsen identification accuracy. Moreover, we must use a precise production calibrated plane. A different perspective approach is to reach EE from points of the same surface, equation, and parameters, of which are known [43]. This surface can be produced by means of 3D printing.
There are also a lot of confidential commercial solutions that are not published in academic journals and proceedings of conferences. These solutions are widely used in conveyers of automotive production. For example, inspection cells perform automated inspection of car bodies with high accuracy, by means of IR. To perform inspection, they have optical sensors with structured light or laser scanners. During work, IR has short and long calibration phases, which can be done in automatic mode. The general idea of such a procedure follows. Calibration artefacts are widely used. These artefacts are constructions with precisely produced spheres, pyramids and so on. Positions of, for example, sphere centers are measured by means of laser trackers. Then, after 10-20 measurements of about 30-150 points of the sphere surface by robot with different sensor orientation, and calculating its center, correction of robot parameters is proceeded by means of least mean squares method or RANSAC. Despite automatic mode, this approach is highly dependent on high-precision location of the equipment, requiring lengthy configuration and tuning, and has high cost of implementation.
Therefore, to eliminate disadvantages outlined above (the need to use precise produced calibrated equipment (plates, probes, spheres, pyramids, etc.), and external high-precision sensors for measuring position of EE in BCS) the similar procedure to the conventional method of calculating tool center point (TCP) can be applied. The data from sensors of rotation angles of actuators in all joints of IR are used in the presented procedure. This data corresponds to the rotation angles of all IR joints when EE IR reaches the same point in space with different orientations, by using a tech pedant. The method of identification of IR kinematic parameters, based on information about IR joints rotation angles will be developed.
This paper presents an essential decrease in procedural cost for identification of IR kinematic parameters, and makes it easily accessible to a wide range of specialists (scientific staff, R&D offices, students, system integrators, production, etc.). This procedure is performed only using two probes and doesn't need expensive external measurement devices, and additional precision equipment. These probes are used for the measurement and generation of an abundance of data for the identification of IR kinematic parameters with quality no worse than that of factory calibration. The relative simplicity of this procedure and absence of hard requirements of equipment allow us to perform identification directly on production lines. Note that the paper constitutes an extension of the conference paper [44], both theoretically and experimentally.

Problem Definition
In the present article, IR with six revolute joints and serial kinematic schemes are considered the most popular type of IR. But the proposed method can be used for any type of IR with a serial kinematic scheme. The kinematic model is described by the following equation [45]: where  k = 1, 6 is the matrix of Denavit-Hartenberg parameters; k is the joint number; Q = (q 1 , . . . , q 6 ) T is the vector of generalized coordinates (rotation angles of joints); is Denavit-Hartenberg matrix [45]. The Denavit-Hartenberg representation forms a matrix of homogeneous transformations (Equation (2)) is 4 × 4 in size, and describes the position and orientation of coordinate system of the k-th joint respectively, the (k − 1)-th joint. For the definition of Equation (2), the following parameters are used: d k is the distance between axis z k−1 and z k along the common normal; θ k is the rotation angle around z k−1 from x k−1 to x k ; a k is the length of common normal; α k is the rotation angle around common normal from z k−1 to z k .
The model (Equation (1)) is used for solving direct kinematic problem, i.e., for calculation of position and orientation of the IR flange in BCS on the base of values, Q and Φ. Moreover, for planning IR movement in BCS, the inverse kinematic problem must be solved, i.e., calculation of desired vector Q * corresponding desired position and orientation T * f of IR flange: where F IKT () is function describing solution of inverse kinematic problem. Equation (3) can be implemented as an analytical expression [46] or a numerical optimization algorithm [47]. Controllers of IR use (Equations (1) and (3)) for calculation of current position and orientation of EE and desired rotation angles of IR actuators. Usually, matrix Φ of nominal geometric parameters, obtained from technical documentation, are used in these expressions. However, precise kinematic parameters Φ of specific IR can be differed from their nominal parameters Φ on some small values because of inaccuracies from fabrication and connection of mechanical parts: Using Φ in Equation (3) leads us to generate vector Q * which differ from vector Q * corresponding desired value T * f on small value β = Q * − Q * . Presence of β leads all actuators to reach positions corresponding Q * and IR flange reaches position T * f which differ from desired position T * f on value ε = T * f − T * f . It means that difference of IR kinematic model parameters from real kinematic parameters of IR leads to errors in position and orientation of IR flange and EE. This is essentially important for the case when trajectory of IR movement is generated on the base of information receiving from external measurement systems (laser or optical scanners, photo cameras, and other) [7][8][9][10][11][12][13][14].
Improvement of IR kinematic parameters can be made by means of special measurement systems, providing high accuracy measurement of linear and angular coordinates of an end effector. But, using such a system is often impossible because of its high cost. At the same time, each IR has a high accuracy measuring system for rotation angles of its actuators that can be used for the calculation of IR kinematic parameters. Thus, the following task is set and solved in this article. The six-degree IR, with serial kinematic scheme, is considered. Its real kinematic parameters are described by matrix Φ. The controller of IR uses matrix of nominal parameters Φ to solve direct and inverse kinematic tasks. This leads to an error ε of EE positioning in BCS. It is necessary to develop the method of identification of IR kinematic parameters based on series of measurements of its generalized coordinates to decrease this error.

Method of Identification of IR Kinematic Parameters
Initial data for matrix Φ estimation is n measurement series of vector Q. Each i-th series includes m i vectors Q i , which are formed as result of EE reaching with different orientations of the same point X i . All measurements are made manually by an operator with the help of an IR teach pendant. The X i position is set by a sharp probe arbitrarily located in the working space. EE movement trajectory is not important as all measurements are made in a stationary position. As a rule, for the convenience of visual inspection by the operator, a sharp probe fixed on an IR flange is used during measurement performance (see Figure 1). the method of identification of IR kinematic parameters based on series of measurements of its generalized coordinates to decrease this error.

Method of Identification of IR Kinematic Parameters
Initial data for matrix Φ estimation is n measurement series of vector . Each i-th series includes vectors , which are formed as result of EE reaching with different orientations of the same point . All measurements are made manually by an operator with the help of an IR teach pendant. The position is set by a sharp probe arbitrarily located in the working space. EE movement trajectory is not important as all measurements are made in a stationary position. As a rule, for the convenience of visual inspection by the operator, a sharp probe fixed on an IR flange is used during measurement performance (see Figure 1). Figure 1. IR with serial kinematic scheme. Following notation is used: q1,…,q6 are rotation angles of IR joints; Oxyz is the base coordinate system; Of xf yf zf is flange coordinate system; is the coordinate vector of the IR flange; , is the coordinate vector of tool center point in BSC; XTCP is the vector of coordinates of the tool center point in Of xf yf zf; is coordinate vector of point.
Thus, as result of measurements, is the following data array: For each vector , = 1, , = (1, ) we can assign vector , of coordinates of tool center point in BCS (see Figure 1) which calculated by means of Equation (1) and matrix : where , ∈ R 3×3 is the orientation matrix of the IR flange in BCS for j-th measurement in i-th series; = 1 , E ∈ R 3×3 is the unity matrix; XTCP is the vector of coordinates of the tool center point in flange coordinate system Of xf yf zf. Thus, as result of measurements, is the following data array: For each vector Q i j , i = 1, n , j = 1, m i we can assign vector X i tool,j of coordinates of tool center point in BCS (see Figure 1) which calculated by means of Equation (1) and matrix Φ: where R i f ,j ∈ R 3×3 is the orientation matrix of the IR flange in BCS for j-th measurement is the unity matrix; X TCP is the vector of coordinates of the tool center point in flange coordinate system O f x f y f z f . Coordinates X i tool,j , calculated by means of Equation (6), differ from coordinates of the real position of the tool center point because of the difference between used IR parameters and their real values. As EE in each measurement series reaches the same point with unknown coordinates, the real values of coordinates of TCP in the same measurement series coincided. This fact can be used for the identification of IR kinematic parameters. Estimation of matrixΦ of IR parameters can be made by adjusting these parameters so coordinatesX i tool,j calculated by Equation (6) move closer to minimal distances. For that reason, we can use following the cost function for this estimation: Equation (7) doesn't include real coordinates of points X i in this case for estimation of IR kinematic parameters, the high-precision measurement equipment is not needed. So, the task of identification IR kinematic parameters is formulated as follows: The method of numerical optimization of Levenberg-Marquardt [48] is used for the estimation of IR kinematic parameters. In this case, the initial data must be presented as follows: Cost function (Equation (7)) can be rewritten as follows: The matrix Φ of IR parameters can be presented as follows: From Equation (11), we can see that 27 parameters are estimated: 24 parameters describe the IR kinematic model, and 3 parameters describe the position of TCP in the The initial estimate of vectorX TCP , we can calculate using the least squares method [49] using array Ξ. Here, we write the expression binding vector X TCP , vector X f position and matrix R f orientation of the IR flange when the probe touches its counterpart in specific points (for example, X 1 = X prob ): Here, we rewrite Equation (12) as follows: where E ∈ R 3×3 is unity matrix. Equation (13) describes the model for the initial estimate of unknown vectorsX TCP and X prob . To use this model, it is necessary to substitute Φ in Equation (1), calculating position and orientation of the IR flange for each measurement and form the following arrays: The estimate of vectorsX TCP andX prob for each measurement series, we can get by means of arrays X f ,i and Ω i and least squares methods: where The initial estimate of vectorX TCP is calculated by means of expression: Obtained estimate ofX TCP and Φ are initial estimate of vectorθ of IR parameters Equation (11), which was used in first iteration of Levenberg-Marquardt method. Identification ofθ by means of this method performs by following expressions: where τ is iteration number; G = ∂P /∂ϑ ∈ R 3L×27 ; E ∈ R 27×27 is the unity diagonal matrix; µ(τ) is a variable defined by speed of tuning; 0 < η < 1 is the step size adaptation parameter. The condition for finishing iteration process (Equation (17) where δ is a small positive constant. The vector P is updated according to Equation (9) on each algorithm iteration considering of improved vectorθ(τ). As result of the work of Equation (17), the estimate of vectorθ is evaluated. This vector provides convergence of pointsX j tool,i in the i-th measurement series to minimal distances. Use of calculated kinematic parameters in robot controller instead its nominal parameters allow to increase accuracy of EE positioning in BCS.

Analysis of IR Parameters Identifiability and Influence of Measurement Errors on Identification Process
Further, we carry out analysis of features of identification IR parameters. Also in this analysis, we neglect indexes corresponding number of measurement series.
Primarily, we write an expression describing dependanceX tool from kinematic parameters of IR, considering that matrix of Denavit-Hartenberg transformation for the k-th join has following form: where N k ∈ R 3×3 is the orientation matrix of the k-th joint in the coordinate system of the (k − 1)-th joint; L k ∈ R 3 is the position vector of the k-th joint in the coordinate system of the (k − 1)-th joint. This dependence with considering of Equation (1) has following view: Considering Equation (17), we can write: where indexes j and s correspond to different measurements in one series. Considering Equation (2), the term (L j 1 − L s 1 ) in Equation (19) has the following view: Parameter d 1 is absent in Equation (20) and, therefore, it is absent in Equations (19) and (7). As a result, this parameter cannot be identified by the offered method. It should be noted that using an incorrect value of d 1 leads to the movement of the position of TCP along the z axis of BCS (see Figure 1) calculated with using this parameter, compared to its real position. In other words, the IR controller will work in coordinate systems which shifted relatively real BCS positions by small values along the z axis. As coordinates of points of the working area model received from the computer vision system, are recomputed to the coordinates system of the controller, then, the presence of the noted shift doesn't influence the accuracy of reaching EE to the desired point in this coordinate system and form of EE trajectory movement. So, impossibility identification of parameter d 1 doesn't influence accuracy in the performance of technological operations.
Further, we consider features of parameter θ 1 identification. Orientation matrix N 1 presents as a multiplication of rotation matrixes on angles θ 1 , q 1 and α 1 : where N θ 1 , N q1 are matrixes of rotation about axis z of the first joint on angles θ 1 and q 1 respectively; N α1 is rotation matrix around axis x of the first joint on angle α 1 . Let us suppose that longitudinal axis of the first link is parallel of z axis of BCS, i.e., α 1 = 0. This situation is typical for most types of IR. So, the variable L 1 considering this supposal and Equations (2) and (21) can be presented as follows: (2)), and X tool , therefore, in form: Considering Equation (22), Equation (19) can be represented as follows: From Equation (23), we see that turning on θ 1 does not affect value X j tool −X s tool as multiplication by the matrix N θ 1 provides only rotation of coordinates of all points ofX tool on angle θ 1 along the z-axis. This rotation will not affect value of the Equation (7), which makes it impractical to tune specified parameter by proposed method.
Next, consider situation when parameter α k = 0, that is, two joints k and k + 1 have axes of rotation located parallel to each other. In this case, orientation matrices of coordinate systems of these joints in BCS will be described as follows: where c k = cos(θ k + q k ), s k = sin(θ k + q k ). The coordinate vectors of position of k and k + 1 joints are determined by the expressions: From Equation (17), it can be seen that partial derivatives of the vector P against d k and d k+1 , taking into account Equations (19) and (24), will have the form: As it can be seen from Equation (26), the partial derivatives for parameters d k and d k+1 for case α k = 0 have the same values. This means that, in this case, it is not possible to separate these parameters, since they will change in the same manner regardless of their actual values. This situation occurs when calculating parameters d 2 and d 3 for PUMA type IR, as well as for a flange with a tool attached to it.
Next, we will consider the effect on the identification process of tool driving errors to the same point in space with different orientations. These errors appear due to physiological characteristics of the human operator, who cannot visually determine approach TCP to a given point, with an error of less than 0.1 mm. In this case, additional errors will be present in the set of measurements, which do not depend on accuracy of kinematic model adjustment. That is, Equation (7) can be rewritten as: where the indices j and s correspond to different measurements in the i-th series; ζ i j is the tool driving errors to a given point during the j-th measurement in the i-th series;X i tool,j are the coordinates of the tool calculated by Equation (6) using vectorQ i j , corresponding to the angles of IR joints rotation when it is accurately driven to a given point in the j-th measurement in the i-th series. That is, vector Q i j can be represented as Q i j =Q i j + E i j , and the occurrence of errors ζ i j are due to the presence of the value E i j . Obviously, in the presence of ζ i j , Equation (27) will not be equal to 0, even with the fine-tuning of the IR parameters, when ||X i tool,j −X i tool,s || = 0. Next, we consider how tool driving errors will affect the process of identification of IR parameters.
Taking E into account, leading to tool driving errors to a given point, and Equation (2), it is possible to write: where e k is corresponding entry of vector E.
It can be seen from Equation (19) that it is possible to reduce the difference between measurements at the presence of the tool driving errors only if the values L k will be reduced. That is, if this will be possible by means of reducing modulus of IR linear parameters (see Equation (28)). Thus, after the initial reduction of Equation (27) due to turning of the majority (except described above, see Equations (20), (23) and (26)) IR parameters, further reduction of its value, which includes ζ i j , is possible by reducing the linear parameters a k and d k . As a result, if the specified process is not limited, then parameters obtained during the identification parameters will no longer correspond to real ones.
To reduce the influence of ζ i j on the process of identification in vector P (see Equation (9)), it is necessary to include only those measurements that satisfy the condition r i p ≥ Y , where Y is an estimate of the tool driven accuracy to a given point. The value of Y is chosen empirically, and its value, in most cases, can be assumed to be 0.1 mm.

Simulation Results
Numerical simulation was carried out to verify the proposed method for identification of IR kinematic parameters. In the simulation, the Mitsubishi RV-2FB robot was considered, with a PUMA kinematic scheme (see Figure 1).
The matrix of nominal parameters of Dennavit-Hartenberg and the matrix of deviations of these parameters from the nominal ones had the form: To test the proposed method, two arrays of data, Ξ 0 and Ξ ζ , containing four series of measurements were generated. They correspond to the approach of the EE X TCP = [10, 20, 109] T with different orientation to points X 1 = [7, −317, 106] T , X 2 = [−366, 10, 106] T , X 3 = [20, 300, 106] T , and X 4 = [200, 200, 106] T . The first array Ξ 0 , corresponded to the case when there were no errors ζ i j , and the second array corresponded to the case when these errors were formed randomly, while their amplitude did not exceed 0.16 mm (see Figure 2).
Next, using the matrix of nominal parameters Φ and arrays Ξ 0 , Ξ ζ according to Equations (14) To test the proposed method, two arrays of data, and , containing four series of measurements were generated. They correspond to the approach of the EE XTCP = [10, 20,  , 106] T . The first array , corresponded to the case when there were no errors , and the second array corresponded to the case when these errors were formed randomly, while their amplitude did not exceed 0.16 mm (see Figure 2).  .25] T . The calculated data was used to identify kinematic parameters by the proposed method. An initial value of J for the array was 728, and the final value after nine iterations of tuning was 2.5 × 10 −4 . For the array , the initial value was 732, and the final value after four iterations of tuning was 8. 16.
As a result of tuning, the following IR parameter values were obtained: As a result of tuning, the following IR parameter values were obtained:  As can be seen from the presented results, without tool driving errors at a given point, the proposed method provides a high accuracy for parameter identification (the error in identifying linear parameters does not exceed 0.01 mm, and angular parameters 0.001 • ). At the same time, if there are specified errors in the measurement array, the accuracy of parameter identification decreases (the error in identification of linear parameters does not exceed 0.07 mm, and angular parameters 0.02 • ).
In addition, from the presented results it can be seen that the parameters d 2 and d 3 have the same value, and their sum corresponds to the real value of the parameter d 2 . Also, the parameters of TCP and the sixth link of the robot are identified together and are not separated. At the same time, the obtained parametersX TCP andφ 6 differ from the reference ones, but they allow us to accurately determine the position of the tool. The z coordinate of the tool, relative to the sixth link, is determined by the sum of d 6 and the z TCP coordinate in the X TCP vector, and the x and y coordinates by the expressions: x = a 6 c 6 + x TCP , y = a 6 s 6 + y TCP . In the real case, when working tools will be installed instead of a probe after identification, it is recommended to carry out a standard procedure for determining TCP, considering the identified parameters. Figure 3 shows the deviation of the TCP coordinates calculated according to the model, Equation (6), using the identified parametersΦ from their true values (solid black line is the deviation when using parametersΦ 0 , gray is when using parametersΦ ζ , and dotted when using nominal parameters Φ).
It can be seen from Figure 3 that identified parameters allow for the calculating position of EE with an accuracy of 0.016 mm, if the measurements were carried out without errors, and with an accuracy of 0.16 mm, if the measurements were carried out with errors, as shown in Figure 2. This accuracy is sufficient to perform most of basic technological operations. Figure 3 shows that when using nominal parameters of IR, accuracy of determining the position of the tool drops sharply, and the error can reach 2.4 mm, which is unacceptable when performing critical technological operations. , = + . In the real case, when working tools will be installed instead of a probe after identification, it is recommended to carry out a standard procedure for determining TCP, considering the identified parameters. Figure 3 shows the deviation of the TCP coordinates calculated according to the model, Equation (6), using the identified parameters from their true values (solid black line is the deviation when using parameters , gray is when using parameters , and dotted when using nominal parameters ). It can be seen from Figure 3 that identified parameters allow for the calculating position of EE with an accuracy of 0.016 mm, if the measurements were carried out without errors, and with an accuracy of 0.16 mm, if the measurements were carried out with errors, as shown in Figure 2. This accuracy is sufficient to perform most of basic technological operations. Figure 3 shows that when using nominal parameters of IR, Thus, the results of the simulation confirmed operability and effectiveness of the proposed method for identification of IR kinematic parameters. It should be noted that use of criterion, Equation (7), makes it possible to obtain 628 differences for the parameter identification from four series of measurements at four stationary points m i = (21,21,18,11). This explains its high accuracy.

Experiment Description
The purpose of the experiment is to verify the method described above. During the verification process, it was necessary to evaluate the accuracy of the calculating position of IR flange center using kinematic parameters obtained during the identification process. For this purpose, a laser tracker was used. It allows us to determine the coordinates of the center of spherically mounted rertroreflectors (SMR) with high accuracy in the associated coordinate system O T x T y T z T . The laboratory setup is shown in Figure 4. Thus, the results of the simulation confirmed operability and effectiveness of the proposed method for identification of IR kinematic parameters. It should be noted that use of criterion, Equation (7), makes it possible to obtain 628 differences for the parameter identification from four series of measurements at four stationary points = (21,21,18,11). This explains its high accuracy.

Experiment Description
The purpose of the experiment is to verify the method described above. During the verification process, it was necessary to evaluate the accuracy of the calculating position of IR flange center using kinematic parameters obtained during the identification process. For this purpose, a laser tracker was used. It allows us to determine the coordinates of the center of spherically mounted rertroreflectors (SMR) with high accuracy in the associated coordinate system . The laboratory setup is shown in Figure 4. When using a tracker to verify this method, there are two main problems. The first is the mismatch of the coordinate systems Oxyz and , which does not allow direct comparison of the coordinates measured by tracker with calculated coordinates of the IR flange position. The second is the mismatch of the center of the SMR, that is, the point whose position tracker measures, with the center of IR flange (see Figure 4), whose coordinates are calculated according to Equation (1) using the identified parameters. These problems were solved as follows.
To perform the comparison of coordinates X FARO measured by tracker with coordinates of flange calculated by model, Equation (1), using parameters identified by the proposed method, it is necessary to make a series of measurements, after which, using the ICP (Iterative Closest Points) algorithm, set of coordinates (cloud of points) should be aligned with cloud X FARO . This will provide an alignment of IR coordinate system with the tracker. Accuracy of the IR kinematic model containing identified parameters will be estimated by the standard deviation between the points of these clouds after their alignment. The smaller this value, the more precisely the kinematic parameters of the IR are estimated.
In order to exclude the influence of mismatch of SMR and flange central point, all movements of the IR during the measurement were carried out with the same flange orientation. In this case, the measured coordinates of SMR will always be shifted by the same value, relative to the center of flange. When using a tracker to verify this method, there are two main problems. The first is the mismatch of the coordinate systems Oxyz and O T x T y T z T , which does not allow direct comparison of the coordinates measured by tracker with calculated coordinates of the IR flange position. The second is the mismatch of the center of the SMR, that is, the point whose position tracker measures, with the center of IR flange (see Figure 4), whose coordinates are calculated according to Equation (1) using the identified parameters.
These problems were solved as follows.
To perform the comparison of coordinates X FARO measured by tracker with coordi-natesX UR of flange calculated by model, Equation (1), using parameters identified by the proposed method, it is necessary to make a series of measurements, after which, using the ICP (Iterative Closest Points) algorithm, set of coordinatesX UR (cloud of points) should be aligned with cloud X FARO . This will provide an alignment of IR coordinate system with the tracker. Accuracy of the IR kinematic model containing identified parameters will be estimated by the standard deviation between the points of these clouds after their alignment. The smaller this value, the more precisely the kinematic parameters of the IR are estimated.
In order to exclude the influence of mismatch of SMR and flange central point, all movements of the IR during the measurement were carried out with the same flange orientation. In this case, the measured coordinates of SMR will always be shifted by the same value, relative to the center of flange.
Thus, experimental studies of the method of identification for IR kinematic model parameters consisted of the following stages.
At the first stage, data was collected. Twenty-one measurements were performed manually at three fixed points in space X 1 -X 3 (see Figure 5) and array Ξ. At the second stage, the kinematic parameters of the considered robot were identified based on obtained data. That is, the matrixΦ was evaluated. Thus, experimental studies of the method of identification for IR kinematic model parameters consisted of the following stages.
At the first stage, data was collected. Twenty-one measurements were performed manually at three fixed points in space X 1 -X 3 (see Figure 5) and array Ξ. At the second stage, the kinematic parameters of the considered robot were identified based on obtained data. That is, the matrix was evaluated. At the third stage, positions of SMR X FARO , mounted on IR flange, were measured in by means of a laser tracker (see Figure 4). Then, the flange was moved, and measurements were performed again. Since SMR was not fixed in the geometric center of flange, its orientation remained unchanged, that is, only linear movements of flange were performed. Simultaneously with recording positions of SMR, the following data was also stored: positions of flange center X UR in BCS Oxyz, position calculated by IR controller, and corresponding values of generalized coordinates QUR. Based on the QUR, according to Equation (1), coordinates of flange center in BCS Oxyz were calculated with the help of parameters from the technical documentation , as well as coordinates evaluated using parameters , identified by the proposed method. Thus, during the measurement process, four point clouds X FARO , X UR , , were obtained. At the fourth stage, for the convenience of further analysis and simplification of the procedure of coordinate system alignment using the ICP algorithm, all geometric centers of point clouds were transferred to the beginning of the coordinate system. To do this, the coordinates of all points in each cloud were recalculated using formula where xp is the coordinates of a specific point in the point cloud; Px is the number of points in this cloud. Finally, at the fifth stage, point clouds were aligned using the well-known ICP algorithm. Moreover, clouds X UR , and were moved to X FARO . Then, by the deviation of points of three clouds X UR , and from points X FARO closest to them and by the total standard deviation, one can indirectly conclude quality of identified parameters. The smaller standard deviation of points after the aligning procedure, the less these points deviate from accurately measured X FARO and the more accurate the kinematic parameters of model, Equation (1), were identified. At the third stage, positions of SMR X FARO , mounted on IR flange, were measured in O T x T y T z T by means of a laser tracker (see Figure 4). Then, the flange was moved, and measurements were performed again. Since SMR was not fixed in the geometric center of flange, its orientation remained unchanged, that is, only linear movements of flange were performed.
Simultaneously with recording positions of SMR, the following data was also stored: positions of flange center X UR in BCS Oxyz, position calculated by IR controller, and corresponding values of generalized coordinates Q UR . Based on the Q UR , according to Equation (1), coordinates X UR of flange center in BCS Oxyz were calculated with the help of parameters from the technical documentation Φ, as well as coordinatesX UR evaluated using parametersΦ, identified by the proposed method. Thus, during the measurement process, four point clouds X FARO , X UR , X UR ,X UR were obtained.
At the fourth stage, for the convenience of further analysis and simplification of the procedure of coordinate system alignment using the ICP algorithm, all geometric centers of point clouds were transferred to the beginning of the coordinate system. To do this, the coordinates of all points in each cloud were recalculated using formula where x p is the coordinates of a specific point in the point cloud; P x is the number of points in this cloud. Finally, at the fifth stage, point clouds were aligned using the well-known ICP algorithm. Moreover, clouds X UR , X UR andX UR were moved to X FARO . Then, by the deviation of points of three clouds X UR , X UR andX UR from points X FARO closest to them and by the total standard deviation, one can indirectly conclude quality of identified parameters. The smaller standard deviation of points after the aligning procedure, the less these points deviate from accurately measured X FARO and the more accurate the kinematic parameters of model, Equation (1), were identified.

Results of Experimental Research
During experimental studies, six degrees-of-freedom IR UR10e with serial kinematic scheme, other than PUMA scheme, and FARO Vantage laser tracker were used. Measurement accuracy of the tracker is up to 0.015 mm (see Figure 4).
In the first stage, data was collected to identify parameters of IR. The probe was driven manually using a teach pendant, and accuracy of this probe driving to given point was inspected visually (see Figure 5). Twenty-one (m 1 = m 2 = m 3 = 21) measurements were performed at three fixed points in space (i = 3). At the same time, vectors Q i j were stored, and thereby, data set Ξ was formed.
At the second stage, identification of IR parameters was performed. The initial value of J, calculated using Equation (7), was 1120, and, after three iterations, it was 28.7. The matrices of identified and nominal parameters had following form:  . Figure 6 shows distances between all pairs of points (there were 630 of them), whose coordinates were calculated based on Ξ using Equation (1), and matrices Φ (gray graph) andΦ (blue graph), respectively. From this Figure, it can be seen that use ofΦ, estimated as proposed in this paper, makes it possible to reduce the distances in each of pairs of points in the X 1 -X 3 measurement series by more than ten times.

Results of Experimental Research
During experimental studies, six degrees-of-freedom IR UR10e with serial kinematic scheme, other than PUMA scheme, and FARO Vantage laser tracker were used. Measurement accuracy of the tracker is up to 0.015 mm (see Figure 4).
In the first stage, data was collected to identify parameters of IR. The probe was driven manually using a teach pendant, and accuracy of this probe driving to given point was inspected visually (see Figure 5). Twenty-one (m1 = m2 = m3 = 21) measurements were performed at three fixed points in space (i = 3). At the same time, vectors were stored, and thereby, data set Ξ was formed.
At the second stage, identification of IR parameters was performed. The initial value of J, calculated using Equation (7), was 1120, and, after three iterations, it was 28.7. The matrices of identified and nominal parameters had following form:  Figure 6 shows distances between all pairs of points (there were 630 of them), whose coordinates were calculated based on Ξ using Equation (1), and matrices (gray graph) and (blue graph), respectively. From this Figure, it can be seen that use of , estimated as proposed in this paper, makes it possible to reduce the distances in each of pairs of points in the X 1 -X 3 measurement series by more than ten times. At the third stage, coordinates of all point clouds (X FARO , X UR , , ) were obtained. Their values are given in Table 1. All of 13 points were measured by the tracker, shown as red dots in Figure 7. For clarity, in the same coordinate system, calculated points are shown: are indicated by blue crosses, are gray circles, X UR are orange plus signs. At the third stage, coordinates of all point clouds (X FARO , X UR , X UR ,X UR ) were obtained. Their values are given in Table 1. All of 13 points were measured by the tracker, shown as red dots in Figure 7. For clarity, in the same coordinate system, calculated points are shown:X UR are indicated by blue crosses, X UR are gray circles, X UR are orange plus signs.  In the fourth stage, geometric centers of all point clouds were calculated, and clouds were moved to the beginning of the coordinate system. The result of this transfer is shown in Figure 8. This figure shows a mismatch of axes of coordinate systems of IR Oxyz and tracker O T x T y T z T , as well as differences caused by use of different sets of IR parameters.
In the fourth stage, geometric centers of all point clouds were calculated, and clouds were moved to the beginning of the coordinate system. The result of this transfer is shown in Figure 8. This figure shows a mismatch of axes of coordinate systems of IR Oxyz and tracker , as well as differences caused by use of different sets of IR parameters. In the last stage, point clouds were aligned by their initial position shown in Figure  8. Its result is shown in Figure 9. This alignment was performed in CloudCompare program, which uses an open PLC library. As a result, standard deviation took the following values: when aligning X FARO and , it was 0.2262 mm, X FARO and X UR was 0.2562 mm, and X FARO and was 0.6715 mm. Deviation of each of the cloud points , X UR , and from corresponding points X FARO is shown in Figure 10. Blue in Figure 10 indicates deviation from X FARO points of cloud , orange is X UR , gray is .   In the last stage, point clouds were aligned by their initial position shown in Figure 8. Its result is shown in Figure 9. This alignment was performed in CloudCompare program, which uses an open PLC library. As a result, standard deviation took the following values: when aligning X FARO andX UR , it was 0.2262 mm, X FARO and X UR was 0.2562 mm, and X FARO and X UR was 0.6715 mm. In the fourth stage, geometric centers of all point clouds were calculated, and clouds were moved to the beginning of the coordinate system. The result of this transfer is shown in Figure 8. This figure shows a mismatch of axes of coordinate systems of IR Oxyz and tracker , as well as differences caused by use of different sets of IR parameters. In the last stage, point clouds were aligned by their initial position shown in Figure  8. Its result is shown in Figure 9. This alignment was performed in CloudCompare program, which uses an open PLC library. As a result, standard deviation took the following values: when aligning X FARO and , it was 0.2262 mm, X FARO and X UR was 0.2562 mm, and X FARO and was 0.6715 mm. Deviation of each of the cloud points , X UR , and from corresponding points X FARO is shown in Figure 10. Blue in Figure 10 indicates deviation from X FARO points of cloud , orange is X UR , gray is .   Deviation of each of the cloud points, X UR , and X UR from corresponding points X FARO is shown in Figure 10. Blue in Figure 10 indicates deviation from X FARO points of cloud X UR , orange is X UR , gray is X UR .
coordinates of points calculated with their help to ones obtained by means of the highprecision laser tracker. That is, the proposed method provided result no worse, and even a little better, than calibration of UR10e robot performed in factory conditions with the help of special precision tools (calibration plate or another calibrated robot). As a result, performed experimental studies have fully confirmed the operability and effectiveness of this method for the identification of IR parameters without use of external measuring devices. Figure 10. Deviation of points from X FARO . Blue indicates deviation from X FARO points of cloud , orange is X UR , gray is .

Conclusions
The paper proposes a method for the identification of kinematic parameters of IR with a series kinematic scheme, which does not require use of external high-precision measuring systems. The proposed method consists of two stages. At the first stage, this is done manually by means of a teach pendant, driven at different orientations to the same fixed point in space and data on the rotation angles of actuators is recorded. At the second stage, using the Levenberg-Marquardt method, the IR model kinematic parameters are tuned in such way as to reduce distances between pairs of tool positions calculated based on mathematical model of this IR. As a result of proposed procedure, it is possible to refine kinematic parameters and, as the research results have shown, significantly increase the accuracy of IR movement in BCS. Experimental studies of the proposed method were carried out and they confirmed its operability and effectiveness.
The proposed method has the following restrictions: it can be used for IR with a serial kinematic scheme, it needs a teach pendant, and an operator, automatization is hardly implemented. However, it can also be used to solve the problem of identifying the elastostatic parameters of IR, which is still relevant. In [50], the results of its use are described, simulation is performed. Experimental research is currently being carried out.  Thus, parameters identified by the proposed method best approximate the coordinates of points calculated with their help to ones obtained by means of the high-precision laser tracker. That is, the proposed method provided result no worse, and even a little better, than calibration of UR10e robot performed in factory conditions with the help of special precision tools (calibration plate or another calibrated robot). As a result, performed experimental studies have fully confirmed the operability and effectiveness of this method for the identification of IR parameters without use of external measuring devices.

Conclusions
The paper proposes a method for the identification of kinematic parameters of IR with a series kinematic scheme, which does not require use of external high-precision measuring systems. The proposed method consists of two stages. At the first stage, this is done manually by means of a teach pendant, driven at different orientations to the same fixed point in space and data on the rotation angles of actuators is recorded. At the second stage, using the Levenberg-Marquardt method, the IR model kinematic parameters are tuned in such way as to reduce distances between pairs of tool positions calculated based on mathematical model of this IR. As a result of proposed procedure, it is possible to refine kinematic parameters and, as the research results have shown, significantly increase the accuracy of IR movement in BCS. Experimental studies of the proposed method were carried out and they confirmed its operability and effectiveness.
The proposed method has the following restrictions: it can be used for IR with a serial kinematic scheme, it needs a teach pendant, and an operator, automatization is hardly implemented. However, it can also be used to solve the problem of identifying the elastostatic parameters of IR, which is still relevant. In [50], the results of its use are described, simulation is performed. Experimental research is currently being carried out. Funding: This research was funded by RFBR as part of the scientific project 20-08-00701, and Sevastopol State University within the framework of an internal grant to attract leading researchers (project identifier 42-01-09/169/2021-5).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
There is no data set associated with the paper.