Joint Dynamics Modeling and Parameter Identification for Space Robot Applications

Long-term mission identification and model validation for in-flight manipulator control system in almost zero gravity with hostile space environment are extremely important for robotic applications. In this paper, a robot joint mathematical model is developed where several nonlinearities have been taken into account. In order to identify all the required system parameters, an integrated identification strategy is derived. This strategy makes use of a robust version of least-squares procedure (LS) for getting the initial conditions and a general nonlinear optimization method (MCS—multilevel coordinate search—algorithm) to estimate the nonlinear parameters. The approach is applied to the intelligent robot joint (IRJ) experiment that was developed at DLR for utilization opportunity on the International Space Station (ISS). The results using real and simulated measurements have shown that the developed algorithm and strategy have remarkable features in identifying all the parameters with good accuracy.


Introduction
Modeling and simulation of the dynamic behavior in a microgravity environment is mandatory for mission success.Not only the reduced gravity is severely changing the dynamic behavior, but also often much more strongly, it is the outer space environment that impacts on physical parameters like joint, structural damping, stiffness in gears and limb structures.That hostile environment is mainly due to big temperature oscillations, solar irradiation, eclipse phases, and hard space radiation.These influences affect predominantly the material, the lubrication properties of space manipulators [1], and other servomechanisms, especially in long-term mission applications.As a result, it is desired a proper knowledge of the time-dependent variances of material behavior in terms of their relevant physical parameters which in turn affects the governing differential equations of motion.Space technology demonstration experiments are required to validate the proposed strategies and algorithms for physical parameter identification.The effect of reduced gravity, temperature [2], and the expected physical parameter change due to material degradation act severely on the proper joint nonlinear dynamic modeling process.These variation especially effects the backlash, friction, stiffness, and control of space manipulators systems.In-flight systems parameters identification, both online and offline versions as well as dynamic model validation are, therefore, a very important pre-requisite to increased confidence in the modeling process.
Slow motion of any mechanical machine has been found to exhibit a highly non-linear friction [3] behavior like: stribeck effect, stick-slip [4], periodic cycle alternating motion, arrest, and so forth.It is well known that the friction and stiffness effect can strongly affect the performance of the robot arm control system, thus, the entire mission success may directly depend on the accuracy of the modeling.In order to obtain a good description of the system, especially in low velocity operation, the nonlinear friction models should be taken into account.However, the identification of nonlinear parameters is extremely difficult to deal with due to the problems of local minima, initial condition, computation time, and so forth.Previous works [3] have reported algorithms that have run time of several days.Such algorithms are almost impracticable if the identification procedure must be performed more than one time, as is the case for space applications, where one is interested in monitoring the parameters time-varying behavior.In this work, a balance between complexity and accuracy is made in order to have a model that accurately describes the friction and stiffness behavior, but also allowing the identification process to be practicable.A friction model that takes into account both, low and high velocity effects, has been derived.The identification strategy uses two versions of LS to identify the parameters, which are linearly dependent upon the measurements.For the nonlinear parameters, a nonlinear global optimization algorithm based on multilevel coordinate search (MCS) [5] has shown a good compromise between accuracy and computation time.

Experiment description
The IRJ experiment (Figure 2.1) developed at DLR-Institute of Robotics and Mechatronics, has served as an experimental setup.The design and construction of IRJ incorporate new features like no bulk wiring on the joint and also a number of sensors that monitor the joint performance.The joints are based on special light-weight harmonic drive (HD) gears, while measuring with high precision all relevant state variables: (a) on the input side, motor angular position and speed via an analogous hall sensor, and commanded electric current, (b) on the output side, off-drive position by using optoelectronics, and a torque measurement device based on strain gauge systems.The sensors used give a high degree of intelligence to the joint.In some tests, two accelerometers have been also attached on the top of the link to measure the acceleration in radial and tangential directions.The motors used are inland brushless DC type, which were redesigned by DLR to provide hollow axes where all cabling are fed through.DLR has also developed a lightweight small robot system with a total weight of less than 20 kg and a length of  1.50 m (Figure 2.2).This design allows a very favorable payload to total weight ratio of about 1 : 3 to almost 1 : 2, compared to conventional industrial robots of more than 1 : 20.
This new design makes the robot very attractive for space-based demonstration missions as on ISS.Currently, there are some studies underway to contemplate about the space experimental use and possible accommodation opportunities at the ISS.However, if not the entire robot system is likely to be operational in the ISS early opportunity utilization phase, the IRJ experiment more probably is expected to get ready for experimental usage.The IRJ experiment will consist of a combination of two of such intelligent rotary joints.The two axes are kinematically combined in order to build up a roll-pitch configuration (Figure 2.2).

Joint dynamics modeling
The main emphasis of the intended space-based identification experiments is directed towards obtaining modeling confidence by proper knowledge of the time-varying joint dynamics parameters, mainly viscous damping, friction/stiction effects, and elasticity within the gears, all of those expecting to be of strongly nonlinear nature.Therefore, the following investigations have been restricted to the modeling and understanding of the nonlinear dynamics of one single intelligent joint.More complex models have already been elaborated for a two-joint configuration and also multibody models have been developed for the seven joint configurations, that is, the entire robotic system, using multibody [1] software code for model generation and simulation.Appropriate identification algorithms have been studied [6,7] and others are still being developed and are underway for these multidegree of freedom systems.
The mathematical model to be used in the identification process is based on Newton's laws that are used to determine the dynamic force interactions and to derive the equations of motions of the joint.Here, only the main steps of the derivation are focused, a detailed description of the modeling is found in [6].
In the joints in the IRJ experiment, the wave generator (wg) is driven by a motor mounted to the circular spline (cs) and the flexspline is attached to the ground.The output is driven by the circular spline.Damping torques, both at the input and output side, have been considered.According to Figures 3.1 and 3.2, while making use only of the pitch (θ) rotary joint, the equations of motion for the IRJ can be described by where T m = K m I a is the applied motor torque with motor constant K m and electric current I a .J in and J out are the input and output inertia, θ in and θ out the respective angular positions.The elasticity within the HD gear is given by the stiffness torque T stiff with T wg = T stiff + T d wg on the input side of the gear and T cs = (N + 1)T wg .N is the gear reduction.
The various damping torques are denoted by T d , attributed with appropriate indices.The applied load on the link (Figure 2.1) side is due to gravity and is given by T load = T g sinθ out with the load amplitude T g .
Based upon the HD manufacturer's catalog values, this gear type typically exhibits the well-known nonlinear behavior.Usually, the dependency between applied torque and the relative angular position Δθ (θ in − θ out ) is given by a combination of piecewise linear functions, T stiff = f (Δθ), depending upon the operational range of the gear.For the identification algorithm to be developed further, it is necessary to replace this piecewise linear behavior by a continuous curve.As a first approach, it has been proved sufficient to apply a third-order polynomial to represent the stiffness torque given by with coefficients k 1 and k 2 to be adjusted.

Wave generator
Gear model Circular spline Moreover, regarding the mechanical nature of the torque measurement system with strain gauges attached to spokes and rings, it may be worthwhile to account also for some compliance in that system.This is necessary to model it as a further spring, being serially connected to the HD gear spring.In total, this would result in a combined softer spring, and can be considered within new stiffness constants k 1 and k 2 that now would enter as unknown parameters within the identification algorithm.
According to experimental results of many authors [3,6], the damping torques T d that appear on the input side, the output side, and inside the HD gear are assumed to capture two facets of damping behavior, namely T visc and T fric .These are a viscous and a dry friction or Coulomb-type part.Thus, total damping torques is written as the viscous part can be strongly nonlinear with a cubic relationship in the angular velocity, with the linearly depending coefficients b 1 and b 2 .For the dry friction, a modified classical Coulomb friction model is required.This is necessary to account for the well-known Stribeck effects.This observes the fact that for low velocities, the friction torque is normally decreasing continuously with increasing velocity, not in a discontinuous manner.Another problem that arise in using the classical Coulomb friction model is the discontinuity at zero velocity.In order to account both problems, Stribeck effects and zero discontinuity, an empirical mathematical model has been adopted, where T N is the normal torque, μ is the friction coefficient, ω S is the Stribeck velocity, i = in,out, δ S is the exponential parameter that is commonly taken either as 0.5, 1 or even 2. Another possibility is to let δ S to be identified by the nonlinear part of the algorithm together with ω 1 and ω 2 .
According to [8], periodic variations in the frictional torque might appear in the HD operation.Thus, we have introduced periodic variations in the frictional torque on the HD output, As this relationship indicates, frictional torque fluctuations of amplitude A cyclic complete one cycle every time the flexspline makes one complete rotation relative to the circular spline.To match this model to experimental observations, a phase shift of γ cyclic is also included.In order to obtain a linear dependency of the two parameters, this relationship Adenilson R. da Silva et al. 7 can be easily transformed to with the new linearly depending parameters A 1 = A cyclic cos γ cyclic and A 2 = A cyclic sinγ cyclic , from where A cyclic and γ cyclic can be recovered.
It is self-evident that not all of the envisaged damping torques given in (3.1) are expected to capture both types, that is, viscous damping and dry friction parts.Where appropriate, only linear viscous damping is considered in order to keep the amount of parameters to be identified at a minimum, as well as the complexity of the joint dynamic model.It has to be kept in mind that the final manipulator configuration consists of seven kinematic degrees of freedom, which otherwise would drive the amount of parameters intensively high.Recalling the given kinematic constraints, the various torques in (3.1) can now be formulated in terms of the input and output positions, θ in and θ out , and their respective velocities ( The damping dependent on the position that appears in T d in is related to Dahl effect [4].
It is necessary to include a position-dependent term also on the input side in order to get good agreement between dynamic model and measured data.

Identification model and strategy
In order to identify the dynamic parameters of the robotic joint, (3.8) have been taken as the dynamic model representation for the identification process.The problem of identifying, especially rigid body dynamics parameters of a robot, has been extensively studied and a vast amount of literature can be found [9][10][11].However, these methods have a common idea: the robot is moved along a selected trajectory while the joint motion and torques are measured.Then, the parameters are offline estimated using a standard offline LS-based technique.In addition, most of these works have used an industrial robot as a test bed.The strategy and algorithm proposed in this paper should guarantee to cope with several requirements, like online procedure, ability to track time-variant parameters, possibility to identify parameters with nonlinear dependency with respect to the measurements in fast way at low-computational cost.
For the algorithm development, (3.8) are rewritten in order to set up a linear combination of the unknown parameters, given by the vector Θ, and the known information, given by the measurement vector φ.The parameters that appear in vector Θ are identified by an RLS with variable forgetting [7] factor and the parameters ω 1 and ω 2 which appear inside the matrix φ are identified by the MCS algorithm.The measured signals are θ in and I a on the input side, θ out , θout , and T out on the output side.The respective velocities θin and θ out and the acceleration signal θin are calculated numerically while regarding filtering techniques to remedy bad measurement signals.
The following specific torque functional relationships have been considered for the dynamic model: ) T d fscs = b fscs 1 sign θout + A 1 sinθ out + A 2 cos θ out , (4.3) where T wg = T n wg = T cs /(N + 1) = T out /(N + 1).Then, the identification model in the linear regression format can be described by where The exponential coefficient δ S has been set to 1. Using the model given by (4.5), a prediction of Y is given by For a given discrete measurement time t k , the predicted error to be minimized in LS sense is Using the singular value decomposition (SVD) approach, the excitation level and linear combination in the information matrix is verified: with U and V being the isometric matrices.If some states are not well excited or there exist some linear combination in the φ matrix, the related singular value σ i will have small Adenilson R. da Silva et al. 9 magnitude, close to machine precision.After testing the matrix φ, the initial condition for the recursive algorithm is obtained by standard batch least squares estimation: and by applying (4.9), one obtains Once the initial conditions are obtained, the recursive identification is carried out by using the algorithm described in [7].

Nonlinear optimization: multilevel coordinate search (MCS) algorithm
Two parameters in (4.1) have nonlinear dependency with respect to measurement data; therefore, they cannot be identified by the RLS approach.There exist several methods that can be used; local minimizer or global one.The local minimizer requires a good starting point and sometimes they deliver only a mathematical solution for the problem.In these cases, the parameters have no longer physical meaning.Most of the global algorithms have very hard computational load, making the identification process almost impracticable.In this paper, the MCS algorithm has been used in combination with RLS approach.
The MCS algorithm has a very interesting combination of local and global search of the minimum.Here, we will point out only the basic ideas of the algorithm, the interested reader is directed to the work of Huyer and Neumaier [5].
Consider the bound-constrained optimization problem with finite or infinite bounds With u and v being n-dimensional vectors with components in R := R ∪ {−∞, ∞} and u i < v i for i = 1,...,n, that is, only points with finite components are regarded of a box [u,v] whereas its bounds can be infinite.If all the bounds are set to infinity, an unconstrained optimization problem is obtained.
The MCS algorithm tries to find the minimizer by splitting the search space into smaller boxes.These boxes contain a distinguished point, the so-called base point, whose function value is known.In splitting the boxes, a nonuniform procedure is used.Parts where low values of the function are expected are carefully examined.In order to speed up the computation procedure, the MCS algorithm combines global search (splitting the boxes with large parts) and local search (splitting the boxes with good function values).This gives a good balance between convergence to the global minimum and computation time.
The nonlinear algorithm requires an index of performance (IP) to be minimized.There exists several ways to define IP criteria.In this work, two criteria have been tested:  (ii) Normal mode: (9) using Θ, start the online algorithm; (10) check N err ; (11) if N err > δ, call MCS algorithm and using the latest measurement window, evaluate the new ΘNL coefficient; else ΘNL is still the minimum (no change in the non-linear parameters); (12) takes the next measurement.Working in this way, the proposed algorithm can track in real time variations in the linear parameters and update the nonlinear parameters only when some corrections are required.The schematic diagram of the integrated algorithm is shown in Figure 6.1.

Results
The proposed strategy and algorithms have been tested in two different situations: first, using only the measured information from IRJ; second, a jump in Θ NL has been simulated in order to check the ability of the algorithm in tracking time variations in Θ NL .Besides, in order to have a normalized system, the data and parameters values of the motor side have been translated to link side, meaning that the gear reduction is 1. (N = 1).

Case 1: using measured data from IRJ.
In this test, the measurements are taken from IRJ with time length of one minute.Figure 7.1, in the upper part shows the motor position and velocity by using a triangular trajectory and on the bottom, the link acceleration.In order to have better resolution, only 20 seconds of measurements are shown.
Using (4.5) as a model and the integrated algorithm, all parameters which appear in Θ have been identified.After 13 seconds collecting data, the matrix φ gets full rank and the starting procedure is completed.The nonlinear parameters are identified by using MCS algorithm and the linear ones are identified by a batch LS. Figure 7.2 shows the convergence process of the non-linear parameters.It can be noted that after few interactions, the convergence criterion has been fulfilled and the nonlinear optimization has been stopped.In the plots, there is a period where all parameters have zero value, this period corresponds to the initialization procedure where there is no online identification.The measurements are collected and an analysis in matrix φ is performed sequentially with the nonlinear and batch estimation.The stiffness coefficients and the nonlinear damping are shown in Figure 7.3.The dashed lines are constant values obtained by an offline procedure using all the data available.It can be noted that all the parameters converge to the offline estimation showing the good convergence and robustness of the RLS algorithm.As expected, the parameter related to the cubic stiffness has low rate of convergence.This fact is early observed in the singular values of the information matrix.The related singular value has the smallest magnitude meaning that this parameter is very difficult to identify.Despite of its small excitation, the cubic stiffness parameter converges to the expected mean value (offline estimation).Finally, Figure 7.6 shows the statistical performance of the identification process.It can be observed that the algorithm has good ability in tracking the reference signal.Most of the errors lies below 5%, the peak of the errors (20%) occurs due to the nature of the trajectory (Figure 7.1) used.In the point where the velocity changes the sign, there exists a peak in the torque and the algorithm cannot predict this high torque immediately.

Case 2: simulation of time-variant parameters.
In order to test the integrated linear and nonlinear identification procedure in case of time-variant parameters, a mixed data set has been used: the angular velocity has been taken from the experiment setup and the friction torque is calculated by setting the values of the parameters in (4.1).The parameter δ s has been set to 1 and the other values used are shown in Table 7.1.By using these values, a simulated friction torque (T d in ) has been obtained to be used in the test of the algorithm.In order to simulate time variant system, the parameters have experimented variations at instant 16 seconds and 32 seconds in the linear and nonlinear ones, respectively.At time 16 seconds, the plant output T d in has been recalculated and a jump of 50% in the linear parameters has been set.Immediately, the RLS algorithm is able to notice the changes in the parameters, and from that it can estimete the new parameters values, as shown in Figure 7.7.At this time, the correction in the linear parameters is sufficient to keep the error smaller than the threshold δ.Therefore, the MCS algorithm has been not activated.The dashed-dot lines represent the parameters values before the jump and the dot lines the values after the jump.Figure 7.7 shows that after the initialization procedure, the parameters identified by the RLS part have fast convergence to the nominal values.The jump in the linear parameters is compensated avoiding the nonlinear optimization.When the nonlinear parameters are changed, the linear ones are affected (peaks in Figure 7.7), but according the convergence in the nonlinear one is reached, the linear parameters approach to the correct values.
At instant t = 32 seconds, the nonlinear parameters have been changed by 20% of their initial values.Then, the norm of error increases and the corrections in the linear parameters are not sufficient to decrease it.Thus, the MCS is activated and the nonlinear parameters are recalculated.When the norm of the error decreases, the nonlinear optimization is stopped and only the fast (RLS) part of the algorithm is running.Figure 7.8 shows the behavior of the nonlinear parameters, it can be noted that the algorithm has fast convergence in both situation: in the initialization and when are recalculated.Due to fast corrections in the linear parameters, change in these parameters does not affect the nonlinear ones.On the other hand, the linear parameters are affected when corrections in the nonlinear ones are required.This occurs because the corrections in the non-linear parameters are not so fast.Due to this feature, the procedure presented here has very lowcomputational load allowing one to track time-variant systems, which contain nonlinear parameters.

Conclusions
In this work, the complete model of the robotic joint has been derived.The obtained model takes into account several nonlinearities; as for the stiffness as well as in the friction model.The typical nonlinear behavior of the friction at low velocity has been taken into account.An integrated (independent linear and nonlinear parts) identification algorithm has been derived and tested by using data from IRJ experiment and also a mixed data to simulate time-variant systems.
The results have shown that strategy presented gives excellent precision at very lowcomputational cost; the integrated algorithm is more than 20 times faster than the completely nonlinear counterpart (if all the parameters is to be identified by MCS algorithm alone).This allows an online identification for almost all of the measurement period, except for a short period, when an update in the nonlinear parameters is necessary; the online identification was not possible.The ability in tracking time-variant parameters has been also tested by using simulated data and the results have shown a fast and accurate response to the variations in both set of parameters: linear and nonlinear ones.Another very important feature of the proposed approach is that there is no necessity of initial guess for all the parameters; they are automatically adjusted by the initialization procedure.It is only required to set the boundary for the nonlinear parameters, even though this is not a requirement but save computation time.

Figures 7 .
Figures 7.4 and 7.5  show the rest of the parameters, which appear in vector Θ.It can be noted that all parameters have stable behavior converging to their expected mean values obtained by full batch identification.Finally, Figure7.6 shows the statistical performance of the identification process.It can be observed that the algorithm has good ability in tracking the reference signal.Most of the errors lies below 5%, the peak of the errors (20%) occurs due to the nature of the trajectory (Figure7.1)used.In the point where the velocity changes the sign, there exists a peak in the torque and the algorithm cannot predict this high torque immediately.

Table 7 .
1. Offline estimation for the parameters.