Aurélie ; Mahieux
email: aurelie.mahieux@ifsttar.fr

Odile Heddebaut

Odile Heddebaut
email: odile.heddebaut@ifsttar.fr

Researcher

Odile Aurélie ; Heddebaut

Aurélie; HEDDEBAUT, Odile Mahieux

Connecting mobility services and spatial territory typology: an application to a former coal mining area in France

Keywords: Mobility services, Land-use, Transport policy, Household Travel Surveys, Transport behaviours, Sustainable mobility. Connecting mobility services and spatial territory typology

When studying "urban unities" as qualified by The French National Institute for Statistics and Economic Studies (INSEE), it may exist different types of towns. They can be centre, periurban or rural areas which could require different transport services. This distinction is based on the percentage of active residents working in an urban centre with an economic approach or upon physical features such as land which has been built and their distances between buildings. Our research field consists in a wide Urban Transport Authority, the "Syndicat Mixte des Transports (SMT) Artois-Gohelle" that counts 115 towns that are globally qualified as urban but they do not share the same degree of transport accessibility. However, they present very different densities of inhabitants.

After defining a new typology of towns the paper questions the possibility to adopt transport policies, improving the whole mobility in a context of sustainable development, to specific travel behaviours by connecting transport modes (TM) in these different areas.

These behaviours, such as commuting to work or travelling for leisure purposes, are obtained by the Household Travel Surveys (HTS) realised in the SMT Artois-Gohelle.

INTRODUCTION

Urban and rural territories are often opposed [START_REF] Bonerandi | Les espaces intermédiaires, forme hybride : ville en campagne, campagne en ville ?[END_REF] and a new category called peri-urban area has emerged (Roux and Vanier, 2008). It raises the issue of transport policy adaptation to each area according to a good knowledge of mobility behaviours. The objective of this paper is to determine a specific typology in order to understand if the use of different transport modes (TM) or innovative services can fit these urban, periurban or rural territories. We distinguish between different TM and allocate them to the different categories of territories. Using a quantitative approach, we analyse the two Household Travel Surveys (HTS) made on our field research. We demonstrate that differentiated mobility behaviours exist among our typology.

ABOUT THE STUDY AREA

Our study area is the Syndicat Mixte des Transports (SMT) Artois-Gohelle. It represents the local authority in charge of transport and mobility policies for its 115 towns. This zone is a former coal-mining area, located in the Nord-Pas-de-Calais Region, in the northern part of France. It has a total of 594,017 inhabitants in 2011 and a surface area of 76,115 hectares that denotes the urban transport perimeter. Few studies exist concerning this area. One study shows a very low household motorization in Lens [START_REF] Lambert | Analyse locale de la motorisation, CREDOC, Collection Des Rapports[END_REF]. It is a polycentric territory with two main centres: Lens with 36,120 inhabitants in 2008 and Béthune with 25,697 inhabitants in 2008. It also presents suburban, peri-urban and rural belts. During the mining era, that ended in 1990, services were concentrated around the mine shaft. Miners and their families could access all the services and jobs on foot [START_REF] Froger | Vieillissement, changement social et pauperisation: le parc de logements face au renouvellement de la population du bassin minier du Nord-Pas-de-Calais (France), Espace Populations[END_REF]. This explained short distances between their homes and workplaces. But now, it is more difficult because jobs, services and shopping centres are more scattered around the territory. Thus, many trips have been created mostly carried out by private car (63% around Lens, 71% around Béthune against 60% for the French average). Nowadays, it is a redeveloping territory carried out by a lot of urban planning projects such as the Louvre-Lens [START_REF] Bodéré | Louvre-Lens: le Louvre autrement, CeROArt. Connecting mobility services and spatial territory typology: an application to a former coal mining area in France MAHIEUX[END_REF], a museum that will be linked with the famous Louvre in Paris to present permanent collections. This aims at giving a new image of this territory. This territory has two important buses with a high level of service (BHLS) projects, the first from Lens-Liévin-Hénin-Carvin and the second from Béthune-Bruay-la-Buissière. They are scheduled on the most frequent bus routes named "BuLLe" with the highest number of passengers. The public transport (PT) network covers all the studied territory with different type of bus services: classical bus for the main cities and demand responsive transport for the more rural towns. In the center of the "SMT Artois-Gohelle", we can imagine that there is a high share for the public transport (because of a high existing supply contrary to the rural zones) and for walking. We can also think that it exists a wide mutual assistance among populations as a legacy of the mining history of this studied zone. A further analysis of the two Household Travel Surveys, available for the territory, allows us to verify or to contradict our assumptions.

DEFINITION OF SPATIAL TERRITORIES

This studied area is mainly considered as an urban one. A review of the literature allows us to consider that it is possible to classify theses territories into urban, peri-urban or rural categories differentiated by their economic [START_REF] Certu | Les modes de transport collectifs urbains, éléments de choix pour une approche globale des systèmes[END_REF], Paquot et al., Brun, 2001, Schmitt and Gofette-Nagot, 2000, Roux and Vanier, 2008), geographic (Poulot, 2008, Bonerand et al., 2004) or sociologic (Thomsin, 2001) approaches and based on density, dwellings or economic activities.

DIFFERENT TYPOLOGIES TO IDENTIFY DIFFERENT SPATIAL TERRITORIES

Our studied area is mainly considered as a dense urban zone. Nevertheless, it is important to obtain different categories of territory to adapt mobility services according to the specificities of each type of area. Four typologies were tested. The first one based on urban unities from INSEE (2010) does not sufficiently distinguish between centre, suburb, isolated and rural towns. Too many towns are in the suburb category. The second typology according to urban areas from INSEE (2010) does not reflect the rural nature of some territories. The third one is a typology according to the population densities. It is not sufficiently discriminating too. Effectively, it only takes into account the habitat without indicating its surface. The last one, described in the below paragraph, is a typology according to the landuse. The aim of this one is to distinguish between the different categories of territories and their classification into the previous definitions explained above. It appears to be the best typology. So, on the following part only this typology is presented. This typology is based on the research of a public land agency in the Nord-Pas-de-Calais Region, called Etablissement Public Foncier (EPF), using habitat data to settle this typology. We add other variables like the share of economic and agricultural activities. We obtain these definitions for eight categories (see Figure 2): -"Centre": that corresponds to the two biggest towns of our study area: Lens and Béthune. Habitat and economic activities represent over 75% of the surface of the towns.

-"Urban pole": that presents the same level of population as the centre, but where economic activities are less prevalent than in the centre (between 23 and 32% against 41% in the centre).

Connecting mobility services and spatial territory typology: an application to a former coal mining area in France MAHIEUX, Aurélie; HEDDEBAUT, Odile 13 th WCTR, July 15-18, 2013 -Rio de Janeiro, Brazil -"Secondary pole": where the level of population is significant compared to our study area, with about one third of the surface of the town devoted to the habitat and a share of economic activities is usually inferior to that of the urban pole.

-"Industrial suburb": where the level of population is inferior to that of the secondary pole but this category counts only a few thousand inhabitants. As the map shows, these towns are located in the former coal-mining area also recognizable by miners'houses due to the former mining activity.

-"Mixed suburb": where the level of population is equal to that of the industrial suburb. The share of agriculture represents more than 50% (between 52 and 71%).

The share of habitats is on average around 25% and the economic activities are still present while agriculture clearly dominates. The mixed suburb towns are mostly located around the centre.

-"Peri-urban": on the edge of our study territory, where the level of population is inferior to that of the industrial and the mixed suburbs. The greatest feature is the dominance of agriculture (between 67 and 99%). Economic activities are almost nonexistent (4% for one town, between 0 and 1% for the others).

-"Mixed peri-urban": that presents the same characteristics as peri-urban. The agriculture still represents more than the half of the town's surface. Economic activities have a stronger weight than that of the peri-urban.

-"Rural": that is located in the south of the study area. The level of population is very low (between 200 and 700 inhabitants, mostly around 300). Agriculture is predominant on these zone's surfaces (between 88 and 94%).

The three first categories are qualified as very high density areas. The two further ones are high density areas and the three later are of low density zones. This typology describes the diversity of the areas within a large urban area. We also apply it to identify specific mobility behaviours, by place of residence, in our study area. We apply this methodology to identify specific mobility behaviour by place of residence.

DIFFERENT TRANSPORT MODES OR NEW MOBILITY SERVICES CROSSED WITH DIFFERENTIATED SPATIAL AREAS

Every TM can be applied to urban territories. Nowadays, several discussions are conducted on the pernicious effects of urban transportation in a context of sustainable development (Bonnafous et al., 1998, Chanei and[START_REF] Chanei | L'environnement dans la décision publique : refonder l'évaluation socio-économique pour des politiques de transport plus durables[END_REF]. The mobility services in rural territories are booming and concrete solutions emerge to facilitate the mobility for people living in these areas (Entreprises Territoires et Développement, 2009). The real issue is to determine which type of mobility services should be implemented in peri-urban areas. The literature also shows what kind of TM fits different distances categories: less than 1 kilometre: walking, 1 to 4 kilometres: bike, 4 to 7 kilometres: two-wheeled motorized, 7 to 10 kilometres: urban PT, more than 10 kilometres: car/ urban PT.

CONNECTING THE TYPOLOGY BASED ON LAND-USE AND TRANSPORT BEHAVIOURS: MAIN RESULTS

Two Household Travel Surveys (HTS) have been conducted on our study area: one for the Lens-Liévin-Hénin Beaumont (LLHB) zone in 2006 and the other for the Béthune-Bruay-la-Buissière-Noeux-les-Mines (BBN) zone in 2005. We have combined these two HTS to sort out the eight different categories of territories described previously. They correspond to the place of residence of people interviewed within the HTS. These territorial categories have been crossed with the distances made by each person and per day and it shows that distances are less important when coming from centre territory (8.8 kilometres). This can be explained by the proximity between dwellings, jobs or services and their concentration in this kind of territory. They increase up to 12.5 kilometres when travelling from mixed peri-urban territories where agricultural activity is present but where industrial activity still exists. Surprisingly, distances are shorter in peri-urban and rural territories where agricultural activity is predominant.

Travel time budget per person and per day is more important where density of urbanisation is much lower. The high level of travel time budget when living in centre territories can be explained probably by traffic congestion effects or by the higher number of trips. In peri-urban and rural areas, travel time is shorter due to this low level of urbanisation avoiding congestion effects and also a lower number of trips.

The main TM used for travelling has been connected within these eight categories of areas. Trips by car increase when the urbanisation degree decreases: for instance, in centre territories, we observe 41% of car use against 73% in rural areas. Secondary pole and industrial suburb present equivalent figures: they have the same share of car driver (respectively 44-45%), of car passenger (20-21%) and of walking trips (respectively 28-26%). The explanation can be an equal proportion of habitat and economic activities in these territories. The difference is only the number of their inhabitants. In the same way, the share of walking trips decreases with the urbanisation degree. A high level of walking when living in centre and urban pole (respectively 35-31%) is observed where the city attractiveness and amenities are more significant, although this level remains high in secondary pole and industrial suburb (respectively 28-26%). Walking trips decline when agricultural activity appears supposing a low level of housings density and greater walking distances such as in mixed suburb and mixed peri-urban territories (respectively 18-14 %).

The predominance of car use in rural territories (73%) is evident, probably for the same reasons. Moreover, walking practice is almost non-existent when living in these peri-urban and rural areas. The share of bike is quite the same (around 2%) in high density zones.

Surprisingly, this share of urban PTs (2%) is very low compared to other French agglomerations such as the nearby urban community of Lille (9%). This share remains constant whatever the place and characteristics of residence territories. The high level of car passenger, particularly in industrial and mixed suburbs (21-22%) compensates also the weak performance of the PT network. It can probably be explained by the mining history and its traditional solidarity.

When we consider the eight categories of territories and the main TM used combined with the travel purpose, we find that car driver has an important share for commuting to work for these eight territories (from 72 to 82%) and for home-accompanying reason (57% in urban poles, 72% in mixed peri-urban areas). Car passenger is very low for commuting to work trips Connecting mobility services and spatial territory typology:

an application to a former coal mining area in France MAHIEUX, Aurélie; HEDDEBAUT, Odile 13 th WCTR, July 15-18, 2013 -Rio de Janeiro, Brazil (9% in peri-urban, 5% in centre). The number of car passengers is higher in centre and urban pole territories for visiting family and friends (25-27%). This motive can effectively be shared by several members of a family. In low-density areas, car passenger rate is higher for hometo-school in peri-urban and rural zones (respectively 42-35%) supposing that walking distances are higher to reach the school or university.

The home-to-school or university trip motive involves walking when living in the highest density zones. More than half of these trips are made on foot (56% in centre and 59% in urban pole). The purposes are similarly represented in industrial and mixed suburbs.

Walking is also involved in leisure activities: in centre (55%) where leisure facilities are great in number and in secondary pole (37%) where they are also present. The rate of walking for leisure activities is high in rural zone (60%) where the environmental and nature amenities favour recreational walks.

Urban PT is mostly used for home-to-school or home-to-university trips (10% in centre areas, 13% in very low-density territories). It can be explained by the age of users who usually do not have a driving licence, by specific services organised by PT for this type of trips and moreover by the fact it is free for school students living in these areas.

Walking is used mainly to go shopping in centres (41%). This rate declines with the degree of urbanisation from 31% in secondary pole and 22% industrial suburb to 6% in periurban zones. As described above, shopping attractiveness is mainly provided in city centres.

For each of the eight territories, we identify the proportion of trips according to their distance from less than 1 to more than 30 kilometres. Everywhere, a great proportion of trips (from 34.1% in the mixed suburb to 53% in the peri-urban) covers distances of less than 1 kilometre. In addition, in centres, urban poles, secondary poles and industrial suburbs twothirds of trips are less than 2 kilometres. Distances travelled between 1 to 4 kilometres represent the second greatest interval: 42% on average for very high density areas, 37% on average for high density areas and 26% on average for low-density areas. Distances travelled between 4 to 7 kilometres represent on average only 8.2% of trips in the very high density areas, 13.3% in the high density areas (also 13.3% in mixed peri-urban territories), and 4.8% in the low-density areas. Using the car more collectively for commuting to work could be proposed for our study area. Effectively, we can deduce from this table that there is some room for enhancing a more collectively use of the car for commuting to work if they are well organised for car passengers. In the same idea, organising walking bus could reduce the use of car to accompany scholars.

For distances under 1 kilometre, walking is the main TM when living in centre and urban pole areas (65-61%). Its share decreases with the urbanisation degree from 65% in centre to 13% in peri-urban and rural territories. This confirms that in high density territories it is easier and pleasant to walk. But these results also show that there is a possibility to enhance walking practice in mixed suburb and mixed peri-urban by offering facilities. Car use for short distances could even be lower in centres, urban poles and secondary poles with an increase of walking amenities for these short distances. When sorting out data scattering TMs by travelled distances, it confirms that for distances of less than 1 kilometre, walking is the dominant TM in the denser zones.

Urban PT is mostly used for distances between 10 to 15 kilometres (10% in centre territories, 8% in urban poles and 11% in mixed peri-urban areas). Despite the existence of PT supply, in centre and urban pole for short distances people rather walk than use the PT possibilities. The bike is mainly used for distances travelled of less than 4 kilometres (of which 1 to 2 kilometres: 3% in centres and secondary poles, 4% in industrial suburbs, 6% in mixed suburbs, 2 to 3 kilometres: 3% in centres and urban poles, 4% in mixed peri-urbans, 3 to 4 kilometres: 5% in peri-urban zones).

The two-wheeled motorised vehicles are also mainly used for distances of less than 4 kilometres (from 2 to 3 kilometres: 3% in centres and mixed suburbs, 4% in peri-urban and 2% in mixed peri-urban areas, and from 3 to 4 kilometres: 3% in urban poles and peri-urban areas and 4% in secondary poles, an exception is made for industrial suburbs where these trips are mainly made from 5 (3%) to 10 kilometres (3%)).

For the same short distance, the second TM is car. In contrast to walking, its share raises when the urbanisation degree decreases (34% in urban pole, 62% in mixed peri-urban, 76% in rural areas). Car use predominates even over short distances. For distances of more than 1 kilometre, cars are principally the distinctly TM. Distances realised by car drivers mainly represent less than 2 kilometres in almost all of places of residence (50% in centres, 52% in industrial suburbs, 55% in peri-urban areas), except for mixed suburbs (36%). For the same interval of distance, car passenger has almost the same share, except for peri-urban territories with 71%.

These results confirm that for long or short distances, trips are mainly made by car.

CONCLUSION

In this paper, a typology of towns is created in order to differentiate territories in a zone globally qualified as urban. This method is based on the distribution of habitat, economic activities and agriculture for each town of our study area.

Technical solutions to improve a more suitable mobility while being less focused on car use were given. However, we do not analyse incentives to implement in order to adopt new mobility patterns. Later, it could be very interesting to focus on the trip chain and to know the origin and the destination of each trip according to the place of residence.

In the full paper we analyse some policy implications for implementing TM in differentiated spatial zones according to our study area. This analysis has shown in particular how industrial restructuring as described in our typology can change both journey lengths and modes.

Source:

  Etude pour l'élaboration du PDU du SMT,SMT, Transétude, 2009.

Figure 1 -

 1 Figure 1 -The "SMT Artois Gohelle" area within Nord-Pas-de-Calais Region

Figure 2 -

 2 Figure 2 -Representation of our study area according to land-use

  Connecting mobility services and spatial territory typology:an application to a former coal mining area in France MAHIEUX, Aurélie; HEDDEBAUT, Odile 13 th WCTR, July 15-18, 2013 -Rio de Janeiro, Brazil

Table I -

 I Differentiated spatial territories for differentiated mobility services

	Urban territories	Peri-urban territories	Rural territories	Not specified territories
	Walking
	Mérenne (2008)
	Walking school bus
	Dupeau (2008)
	Mackett et al.

ACKNOWLEDGEMENTS

This research has been funded by the Nord-Pas-de-Calais Region and the French Environment and Energy Management Agency (ADEME) within the research project SUIM (Innovative Urban Services for a new Mobility). The authors wish to acknowledge these two institutions for their helpful support. We gratefully acknowledged the detailed comments of two anonymous referees and WCTR 2013 conference participants for their helpful advices.


A Klimchik

Y Wu

S Caro

C Dumas

B Furet

A Pashkevich
email: anatol.pashkevich@mines-nantes.fr

Modelling of the gravity compensators in robotic manufacturing cells

Keywords: industrial robot, calibration, gravity compensator, geometrical modeling. 

The paper deals with the modeling and identification of the gravity compensators used in heavy industrial robots. The main attention is paid to the geometrical parameters identification and calibration accuracy. To reduce impact of the measurement errors, the design of calibration experiments is used. The advantages of the developed technique are illustrated by experimental results .

INTRODUCTION

Currently, aeronautic industry requires high-precision machining of huge aircraft component made from the high performance materials. For these applications, robots are quite attractive due to their large workspace that can be also easily extended. Another considerable advantage is their capability to process the parts with complex shape and geometry. However, processing of these materials causes essential cutting forces that are generated by tool-workpiece interaction during the material removal. Since these forces decrease the machining accuracy, the robot manufactures try to improve the robot stiffness by increasing the link crosssections. This approach obviously leads to augmentation of robot link weights. So, the gravity forces applied to the manipulator components become non-negligible and also contribute to the position errors. To overcome this difficulty, the link weights are tended to be balanced by gravity compensators, which considerably complicate the stiffness modeling of these heavy manipulators.

The problem of stiffness modeling for the heavy manipulators with gravity compensators has been in the focus of rather limited number of works. In contrast, for conventional serial manipulators without gravity compensators, the problem has been studied by a number of authors that considered both industrial and medical robots with essential compliance in the links and joints (Meggiolaro 2005, Kövecses 2007). Relevant works are mainly based on the virtual joint method (VJM), which lumps the elastostatic properties of the robot components in virtual springs. To our knowledge, the stiffness modeling for the manipulators with gravity compensators has not been studied in detail yet. Currently, the main activity in this area focuses on the gravity compensator design (Takesue 2011[START_REF] Luca | A PD-type Regulator with Exact Gravity Cancellation for Robots with Flexible Joints[END_REF]. On the other hand, since the considered robots include closed loops induced by the compensators, some technique previously developed for the parallel manipulators can be adopted [START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot with uncoupled kinematics[END_REF][START_REF] Company | A Method for Modeling Analytical Stiffness of a Lower Mobility Parallel Manipulator[END_REF], Pashkevich 2011). This paper focuses on the geometrical and stiffness modeling of the spring-based gravity compensators that can be integrated in a VJM-based stiffness model of a serial manipulator (Klimchik 2012). The main attention is paid to the identification of the model parameters and calibration experiment planning. The developed approach is confirmed by the experimental results that deal with the industrial robot employed in manufacturing of large-dimensional aircraft components.

To address these problems, the remainder of the paper is organized as follows. Section 2 presents the gravity compensator model. In Section 3, calibration methodology is presented. Section 4 proposes identification algorithms. Sections 5 is devoted to the experiment design. Experimental validation is presented in Section 6. Finally, Section 7 summarizes the main contributions.

GRAVITY COMPENSATOR MODEL

The mechanical structure of the gravity compensator under study is presented in Fig. 1. The compensator incorporates a passive spring attached to the first and second links, which creates a closed loop that generates the torque applied to the second joint of the manipulator. Corresponding model is presented in Fig. 2, where the most essential geometrical parameters are denoted as x a , y , and the compensator is described by the spring compliance c k and the preloadi g 0 a L n s . This design allows us to limit the stiffness model modification by incorporating in it the compensator torque and adjusting the virtual joint stiffness matrix that here depends on the second joint variable only.

q

The compensator geometrical model includes three node points P 0 , P .

Besides, let us introduce the angles  ,  and the distances This mechanical design allows to balance the manipulator weight for any given configuration by adjusting the compensator spring preloading. It can be taken into account by introducing the zero-value of the compensator length 0 s corresponding to the unloaded spring. Under this assumption, the compensator force applied to the node P 1 can be expressed as follows 0

•( ) sc s FK s   (2)
where is the compensator spring compliance.

c k
Further, the angle  between the compensator links P 0 P 1 and P 1 P 2 (see Fig. 2) can be found from the expression

•(1 )• • • ( /s i n cc ss L Ma K    2 ) q 2 q (4)
Upon differentiation of the latter expression with respect to 2 , the equivalent stiffness of the second joint (comprising both the manipulator and compensator stiffnesses) can be expressed as:

q 22 0 •• c KKK L a    (5)
where the coefficient

2 2 0 22 2 sin cos co •( s )( ) ( q s aL qq s s            2 ) q  (6)
highly depends on the value of joint variable 2 and the initial preloading in the compensator spring described by q 0 s .

Hence, using expression (5), it is possible to extend the classical stiffness model of the serial manipulator

  1 (F) 1 (F)T C θθ θ θθ ()    KJ K HJ (7)
by modifying the virtual spring parameters in accordance with the compensator properties. Here, C θ are stiffness matrices in the Cartesian and joint spaces respectively, θ are corresponding Jacobian and Hessian matrices (for more details see Klimchik 2012). While in the paper this approach has been used for the particular compensator type (spring-based, acting on the second joint), the similar idea can be evidently applied to other compensator types.

, KK (F) θθ , JH Summarizing this Section, it is worth mentioning that the geometrical and elastostatic models of a heavy manipulator with a gravity compensator should include some additional parameters ( x a , y , and c a L K , 0 s for the presented case) that are usually not included in datasheets. For this reason, the following Sections focus on the identification of the extended set of manipulator parameters.

MODEL CALIBRATION METHODOLOGY

In contrast to the serial manipulator that can be treated as a principal mechanism of the considered robots, geometrical data concerning gravity compensators are usually not included in the technical documentation provided by the robot manufacturers. For this reason, this Section focuses on the calibration methodology of the geometrical parameters for the described above compensator mechanism (see Fig. 1).

The geometrical structure of the considered gravity compensator is presented in Fig. 2. Its principal geometrical parameters are denoted as , L

x a , y , where

a •cos x aa   , •sin y aa  
(see notation in Section III). As follows from the figure, the identification problem can be reduced to the determination of relative locations of points P 0 and P 1 with respect to P 2 .

It is assumed that the measurement data are provided by the laser tracker whose "world" coordinate system is located at the intersection of the first and second actuated manipulator joints. The axes Y, Z of this system are aligned with the axes of joints #1 and #2 respectively, while the axis X is directed to ensure right-handed orthogonal basis. To obtain required data, there are several markers attached to the compensator mechanism (see Fig. 3). The first one is located at point P 1 , which is easily accessible and perfectly visible (the center of the compensator axis P 1 is exactly ticked on the fixing element). In contrast, for the point P 0 , it is not possible to locate the marker precisely. For this reason, several markers P 0i are used that are shifted with respect to P 0 , but located on the rigid component of the compensator mechanism (these markers are rotating around P 0 while the joint coordinate 2 is actuated). It should be noted that for the adopted compensator geometrical model (which is in fact a planar one), the marker location relative to the plane XY is not significant, since the identification algorithm presented in the following sub-section will ignore Z-coordinate.

q Using this setting, the identification problem is solved in two steps. The first step is devoted to the identification of the relative location of points P 1 and P 2 . Here, for different values of the manipulator joint coordinates 2

{ i , q 1, } im  ,
the laser tracker provides the set of the vectors   1 i p describing the points that are located in an arc of the circle. After matching these points with a circle, one can obtain the desired value of (circle radius) and the Cartesian coordinates 2 of the point P 2 (circle center) with respect to the laser tracker coordinate system. The second step deals with the identification of the relative location of points P 0 and P 2 . Relevant information is extracted from two data sets   01 i p and that are provided by the laser tracker while targeting at the markers P 01 and P 02 . Here, the points are matched to two circle arcs with the same center (explicitly assuming that the compensator model is planar), which yields the Cartesian coordinates 0 of the point P 0 , also with respect to the laser tracker coordinate system. Finally, the desired values   As follows from the presented methodology, a key numerical problem in the presented approach is the matching of the experimental points with a circle arc. It looks like a classical problem, however, there is a particularity here caused by availability of additional data   2i q describing relative locations of the points  . This feature allows us to reformulate the identification problem and to achieve higher accuracy compared with the traditional approach.

 1 i p 4. IDENTIFICATION ALGORITHMS
The above presented methodology requires solution of two identification problems. The first one aims at approximating of a given set of points (with additional arc angle argument) with an arc circle, which provides the circle center and the circle radius. The second problem deals with an approximation of several sets of points by corresponding number of circle arcs with the centers on the same rotation axis. Let us consider them sequentially.

Algorithm #1

To match the given set of points  with additional set of is the orthogonal rotation matrix, is the vector of the translation that defines the circle center. It worth mentioning that such a formulation has an advantage (in the sense of accuracy) comparing to a traditional circle approximation and it is a generalization of Procrustes problem known from the matrix analysis.

R t

Using equation ( 8), the identification can be reduced to the following optimization problem

    1 ,, min
T ii i ii m F           Rt pR u t pR u t (9)
which should be solved subject to the orthogonality constraint T  RR I. After differentiation with respect to t , the latter variable can be expressed as

1 11 1 i m i mm      tp R i m i   u (10)
That leads to the simplification of (9) to

    1 , min m i r T ii ii F       R pR upR u     (11)
where 11 11

;

ii i i m i i i i m m      pp p uu u m   (12)
Further, differentiation with respect to  yields to

11 m T ii i i m i     pR u uu T i    (13)
So, finally, after relevant substitutions the objective function can be presented as

11 1 12 min mm m ii i TT T ii ii i i F                 R pp uu pR u     ( 14
)
where the unknown matrix must satisfy the orthogonality constraint and can be solved using SVD-decomposition of the matrix

1 m i i T i    up UΣ V   ( 16
)
where the matrices are orthogonal and is the diagonal matrix of the singular values. Further, using the same approach as for the Procrustes problem, it can be proved that the desired rotation matrix can be computed as

, UV Σ T  RV U (17)
which sequentially allows to find the scaling factor  defining the arc radius and the vector defining the arc center.

Algorithm #2

The second problem aims at approximating of several point sets i by corresponding number of concentric circle arcs with the centers 0 on the rotation axis n . It should be noted that the data set   1 ,..., k i pp j p   i q is not useful here, since the required angles  i  are not measured directly and cannot be computed without having exact compensator geometry. In this case, the objective function can be written in the straightforward way

      0 2 2 00 11 , min j j km jj jj ji i ji R T FR      p pp pp  (18)
where j R denotes the j-th ark radius and 0 j p is the corresponding center point. However, for this formulation it can be easily proved that the optimization problem (18) does not lead to a unique solution. In fact, it gives the rotation axis passing through the center points 0 j p , which can be expressed as

0 c j j   pp n ( 19
)
where is the axis direction vector, c p is a point belonging to the axis, and n j  are corresponding scalar factors.

To solve the problem (18), first the objective function F can be differentiate with respect to 2 j R that yields the following expressions for the arc radii

   21 0 1 m jj j T  0 j ji i i Rm     pp pp  m i n s  (20)
Further, after relevant substitution, the objective can be rewritten as that after substitution in it ( 19), ( 24) and ( 26) gives

   
  2 11 min km j i j T i f    n pn   (28)
Further, differentiation (28) with respect to , the optimization problem reduces to the solving following homogeneous linear equation system

n 11 km jj ii ji T    pp n 0   (29)
Non-trivial solution of this system can be found using the singular value decomposition of the matrix 11

km j j ii ji T   pp   11 km jj ii ji T    pp U Σ V T   (30)
where the vector is the last column of the matrix . n V

It should be mentioned that practical application of the latter expressions is essentially simplified by the adopted assumption concerning orientation of the reference coordinate system, where the direction of the rotation axis is close to Z-direction.

n Hence, the developed algorithms allows us to identify the compensator geometrical parameters , L x a , y that are directly related to the above mentioned rotation center points P 0 , P 2 and corresponding radii. Below they will be applied to the processing of the experimental data. a

DESIGN OF CALIBRATION EXPERIMENTS

The main idea of the calibration experiment design is a set of robot configurations (as well as marker locations) that ensure the best identification accuracy. The key issue here is the ranging of different plans in accordance with the prescribed performance measure.

For the considered identification problem, the design variables are the set of angles  and the marker locations. The objective functions to be minimized are computed via the covariance matrix that describes the identification errors for the geometrical parameters and to be estimated. Since two identification algorithms are independent, selection of the optimal configurations  2i q L a   2i q and marker locations can be considered sequentially.

Assuming that each experiment includes the additive measurement errors in the Cartesian coordinates i ε , expression (13) allows us to present the variance of the parameter  in the following way

     2 11 1 var TT T T ii mm m ii i i ii i E       uRεε Ru u u    (31)
where denotes the expectation and the orthogonal matrix defines the orientation of the base coordinate system. Following usual assumption concerning the measurement errors (independent identically distributed, with zero expectation and standard deviation var

T ii m i     uu   (32)
Further, it can be proved that

1 i  m T ii m   uu   q L
. So, the variance (32) does not depend on the angles 2i . Thus, the identification accuracy for the parameter depends on the number of experiments only.

For the remaining geometrical parameter a , the identification error depends on the estimation precision of relative location of the points P 2 and P 0 . Since relevant identification algorithms employ independent measurement data, the variance can be computed as the sum of the traces of 2 cov and 0 cov , where 1 t and 2 are the vectors of Cartesian coordinates for the points P 2 and P 0 .

 var a ( ) t ( ) t t
For the point P 2 , expression (10) leads to the following covariance matrix

  11 1 1 22 2 cov m T i mm m ii i i ii i mE         t εε Ru u R  T T  u is simplification is based on the (33)
which can be further simplified down to

  12 2 1 2 1 cov T i mm i ii mm       tI u (34)
Th above derived expressions

  1 2 1 T ii mm ii Em     εε I and
and on the assumption that z-axis of the coordinate system is directed along the second joint axis. Hence, for the point P 2 , the optimization problem that is related to the design of calibration experiment, can be formulated as

  22 / E   m q      2 22 2 1 2 1 cos sin min i i m i i m i q q F     (35)
This problem should be solved taking into account joint limits of the industrial robot. In the case when the range of angles is over 2 q  , it is possible to achieve zero value of this objective since equations are solvable. It should be noted similar equations arise in calibration experiment design for some robots without gravity compensators and have been studied in details in our previous work (Klimchik 2011).



For the point P 0 , similar expression includes a set of the angles i  that can be recomputed to the joint angles requires for the manipulator control (see Fig. 3). Here, it is reasonable to find optimal marker locations on the rigid part of the gravity compensator. It can be proven that using these assumptions, the design of experiment reduces to the following optimization problem 2i q    

EXPERIMENTAL RESULTS

To demonstrate efficiency of the developed technique, the experimental study has been carried out. The experimental setup employed the robot KR-270 and the Leica laser tracker, which allowed us to measure the Cartesian coordinates of the markers attached to the compensator elements with the accuracy of 10 μm (see Figs 1,3). Six different manipulator configurations where considered that differed in the value of the joint angle 2 and three markers has been used. The experimental data are presented in Table 1.

q These data has been processed using the developed identification algorithm presenting in the Section 4. The obtained values for the parameters of interest , L x a , y are given in Table 2, which also includes the identification errors computed using the Gibbs sampling technique. a In addition, there were evaluated the elastostatic properties of the gravity compensator. Corresponding curves describing influence of the compensator on the equivalent stiffness of the manipulator joint (see eq. ( 6)), are presented in Fig. 4. They demonstrate essential non-linearity of the compensator impact throughout of the robot workspace, which, in addition, highly depend on the spring preloading 0 s . -31.84 183.86 -872.10 -125.38 -813.50 -255.59 -30 -118.44 143.42 -872.30 -126.07 -813.33 -256.18 -60 -173.30 65.12 -872.50 -109.90 -825.09 -244.64 -90 -181.76 -30.14 -868.43 -78.20 -844.66 -219.04 -120 -141.45 -116.82 -858.90 -47.60 -859.43 -190.44 -145 -78.10 -165.47 -852.53 -33.68 -864.66 -176.01

0[m] s  0 0.2 [m] s  0 0.2 [m] s  0 0.4 [m] s  0 0.6 [m] s  0 0.8 [m] s  2 ,[deg] q 
The paper presents a new approach for modeling and calibration of heavy industrial robots with gravity compensators. It proposes a methodology and data processing algorithms for the identification of the gravity compensator geometrical parameters. To increase the identification accuracy, the design of experiments has been used aimed at proper selection of the measurement configurations and marker point locations. The advantages of the developed techniques are illustrated by experimental study of the industrial robot Kuka KR-270, for which the model parameters of the gravity compensator have been identified. Another set of experiments have been carried out to identify the elastostatic properties of the considered manipulator. There were considered 15 different configurations, which have been found taking into account physical constraints that are related to the joint limits, work-cell obstacles and safety reasons. To ensure identification accuracy for each configuration, the experiments were repeated three times. In total, 405 equations were considered for the identification, from which 7 physical parameters have been obtained. Corresponding values of the elastostatic parameters for the gravity compensator and for the manipulator are presented in Table 3, where k i denotes the i-th joint compliance. There were also computed the identification errors (using the Gibbs sampling technique). As follows from Fig. 5, the gravity compensator essentially reduces the equivalent compliance of the second joint (compared to the serial manipulator without the gravity compensator). Geometric and elastic error calibration of a high accuracy patient positioning system. Mechanism and Machine Theory, vol. 40, pp. 415-427. Pashkevich A., Klimchik A., Chablat D. (2011). Enhanced stiffness modeling of manipulators with passive joints. Mechanism and Machine Theory, vol. 46, pp. 662-679. Takesue N., Ikematsu T., Murayama H. and Fujimoto H. (2011). Design and Prototype of Variable Gravity Compensation Mechanism (VGCM). Journal of Robotics and Mechatronics,vol. 23,no. 2,. Compliance of equivalent non-linear spring in the second joint

aFig. 2 .

 2 Fig. 1. Gravity compensator of robot KUKA KR-270 TM

Fig. 3 .

 3 Fig. 3. Geometrical parameters of the gravity compensator and location of the measurement points labelled with markers

  of reference points located on the unit circle whose distribution on the arc is similar to i , [cos , sin , 0] T ii i qq  u p  is the scaling factor that defines the desired circle radius,

  Since the matrix is included in the second term only, the problem can be further simplified to R

  t



  are the angles around the point P 0 between the compensator spring and j-th marker location. It is clear that in this case the best solution is produced by similar to (35), this problem can be easily solved by locating the markers on the opposite sides of the compensator rotation axis. reduced to the solution of the above presented trigonometric equations that allows us essentially increase the calibration accuracy.

Fig. 4 .

 4 Fig. 4. Variation of the gravity compensator impact on the equivalent stiffness of the second joint

  , Angeles J.,(2007). The stiffness matrix in elastically articulated rigid-body systems. Multibody SystemDynamics, pp. 169-184. Meggiolaro M.A., Dubowsky S., Mavroidis C. (2005).

Table 1 . EXPERIMENTAL DATA

 1

	q 2	P 1		P 01		P 02
	[deg]	x, [mm]	y, [mm]	x, [mm]	y, [mm]	x, [mm]	y, [mm]
	-0.01

Table 2 . GEOMETRICAL PARAMETERS

 2

		L, [mm]	a x , [mm]	a y , [mm]
	value	184.72	685.93	120.30
	accuracy	±0.06	±0.70	±0.69

Table 3 . ELASTOSTATIC PARAMETERS

 3 Klimchik A., Bondarenko D., Pashkevich A., Briot S., Furet, F. (2012). Compensation of tool deflection in roboticbased Milling. The 9th International Conference on Informatics in Control, Automation and Robotics (ICINCO 2012), pp. 113-122.

	Parameter		value			accuracy
	k c , [rad×μm/N]	0.144		±0.031 (21.5%)
	s 0 , [mm]	458			±27 (5.9%)
	k 2 , [rad×μm/N]	0.302		±0.004 (1.3%)
	k 3 , [rad×μm/N]	0.406		±0.008 (2.0%)
	k 4 , [rad×μm/N]	3.002		±0.115 (3.8%)
	k 5 , [rad×μm/N]	3.303		±0.162 (4.9%)
	k 6 , [rad×μm/N]	2.365		±0.095 (4.0%)
	2.95
	2.9
	2.85
	2.8
	-150 2.75	-120	-90	-60	-30	0


Y Wu

A Klimchik
email: klimchik@mines-nantes.fr

A Pashkevich
email: pashkevich@mines-nantes.fr

S Caro

B Furet

Efficiency Improvement of Measurement Pose Selection Techniques in Robot Calibration

Keywords: robot calibration, measurement pose selection, test-pose approach 

The paper deals with the design of experiments for manipulator geometric and elastostatic calibration based on the test-pose approach. The main attention is paid to the efficiency improvement of numerical techniques employed in the selection of optimal measurement poses for calibration experiments. The advantages of the developed technique are illustrated by simulation examples that deal with the geometric calibration of the industrial robot of serial architecture.

INTRODUCTION

In the usual engineering practice, the accuracy of a manipulator depends on a number of factors. Usually in robotics, the geometric and elastostatic errors are the most significant ones. Their influence on the robot positioning accuracy highly depends on the manipulator configuration and essentially differs throughout the workspace. To achieve good accuracy in all working points, adequate geometric and stiffness models are required. While the model structure is usually well known, the identification of the model parameters (calibration) is rather time consuming and requires essential experimental work. For this reason, optimal selection of measurement poses for robot calibration is an important problem, which is still in the focus of numerous research papers [START_REF] Niversity Press | choosing measurement poses for robot calibration with th convergence method and tabu search[END_REF], Sun 2008).

At present, the main activity in this area is concentrated around the geometric calibration (Khalil 2002). On the other hand, the elastostatic calibration which is also very important for many applications (such as precise machining) has attracted less attention of the researchers (Meggiolaro 2005). However, for both of these calibration procedures, the problem of measurement pose selection is one of the key issues allowing to reduce essentially the measurement error impact [START_REF] Klimchik | Enhanced stiffness modeling of manipulators Mechanism and Machine Theory[END_REF]. At first sight, this problem can be solved using well known results from the classical design of experiments theory. However, because of the specificity and nonlinearity of the manipulator geometric and elastostatic models, the problem solution is not so obvious. The main difficulties here are in the area of definition of a reasonable optimality criterion (which has clear engineering sense) and also in efficient solution of the relevant optimization problem, which has rather high dimension.

Among related works, it is worth mentioning several papers. The majority of the measurement pose selection techniques relies on the optimization of some functions depending on the singular values of the identification Jacobian. For example, Zhuang used genetic algorithm for minimization of the condition number of this matrix [START_REF] Zhuang | Optimal planning of robot calibration experiments by genetic algorithms[END_REF]. In other work [START_REF] Niversity Press | choosing measurement poses for robot calibration with th convergence method and tabu search[END_REF], to decrease the sensitivity to local minima, Daney developed the local convergence method and Tabu search technique based on the observability index. However, the performance measures used in these works are rather abstract and are not directly related to the robot accuracy. Besides, the related objective functions are very difficult for the optimization due to existence of a number of local minima.

To find the global one, heuristic search is usually used as the numerical algorithms, which often require tedious computations. All these motivate the research direction of this work.

In this paper, the problem of optimal design of calibration experiments is studied for the case of robot manipulator of serial architecture. In contrast to other works, the optimization problem related to measurement pose selection is formulated using the proposed performance measure (test-pose approach), which has clear physical meaning and is directly related to robot accuracy. The main attention is paid to the efficiency improvement of the related numerical routines.

PROBLEM STATEMENT

Geometric and Elastostatic Models of Manipulator

In industrial robot controllers, the end-effector position of the manipulator is usually computed using the geometric model. For some specific applications, such as high-speed machining that generate essential external loading, the elastostatic model should be also used. However, in practice, the robot geometric parameters essentially differ from the nominal values declared in technical specifications and vary from one robot to another. In addition, elastostatic parameters of the manipulator are not provided by the robot manufacturers and can be identified from the experiments only. So, the manipulator model parameter identification is an important step in practical application of industrial robots.

The manipulator geometric model provides the position/orientation of robot end-effector as a function of the joint variables and its inherent parameters. This model is usually presented as a product of homogeneous transformation matrices, which after some transformations can be presented as the vector function

 , g  pq П
(1)

where vector denotes the end-effector position, vector aggregates all joint angles and П are the vector of unknown parameters to be identified. These unknowns differ with the applied parameterization methods in robot geometric modelling, such as the classical Denavit and Hartenberg approach and its modified version (Khalil 1986). In this paper, there are considered the most essential components of the vector , which are the deviations of the robot link lengths i p q П l  and the offsets i  in the actuated joints. Since the deviations of geometrical parameters  are usually relatively small, calibration usually relies on the linearized model

q П     00 ,, g g  pq П Jq ПП  /   p   (2)
which includes the conventional geometric Jacobian computed for the nominal geometric parameters .

    00 ,, g g   q П q П J Π 0 П
The elastostatic properties of a serial robotic manipulator represent its resistance to deformations caused by external forces/torques and are usually described by the Cartesian stiffness matrix , which is computed as

C K 1  C θθ θ T   KJ K J (3)
where θ is a diagonal matrix that aggregates the joint stiffness values (that are the unknowns to be identified) and θ is the corresponding elastostatic Jacobian. This model can be derived using the virtual joint method, which describes all elastostatic properties of compliant elements by localized virtual springs located in the actuated joints (Salisbury 1980).

Using the Cartesian stiffness matrix, the elastostatic model (or force-deflection relation) can be expressed as

K J 1 θθ θ T   JJ wK (4)
where is the position deflection at the robot end-effector caused by the external wrench w , which integrates both the external force and torque. This linear relation can be further used for the calibration where the desired parameters to be identified are the components of matrix .

p θ K
In the frame of this work, several assumptions concerning calibration of these models are accepted:

A1: For the geometric calibration, each calibration experiment produces two vectors  ii , which define the robot end-effector displacements and corresponding joint angles.

, pq

The linear relation between the errors in geometric parameters and the end-effector position deviations can be written as = ( )

ig i  pJ qΠ (5)
where ()

g i
Jq is the Jacobian matrix that depends on manipulator configuration i q and vector collects the unknown parameters to be identified.

Π

A2: For the elastostatic calibration, each calibration experiment produces three vectors

, where defines the applied forces and torques.

  ,,

ii i pqw i w

In accordance with (Pashkevich 2011), the corresponding mapping from the external wrench space to the end-effector deflection space can be expressed as

θθ θ = () () T ii i pJ q k J q w i (6)
where θ is a matrix that aggregates the unknown compliance parameters

k   n k

Identification of the Model Parameters

The problem of parameter identification of the robot manipulator can be treated as the best fitting of the experimental data by corresponding models. These data are measured under several assumptions concerning the measurement equipment:

A3: The calibration relies on the measurements of the endeffector position only (Cartesian coordinates   ,,

x yz ppp ).

A4:

The measurements errors i accommodated in each measurement of end-effector position are treated as independent identically distributed random values with zero expectation and standard deviation ε  .

For computational convenience and taking into account the influence of measurement errors, the geometric and elastostatic models described by separate linear equations ( 5) and ( 6) can be expressed in the following integrated form ()

ii i     pB q X ε (7)
              J0 B0 A JA (8)
where is the Jacobian matrix that can be obtained by differentiating the manipulator geometric model with respect to the desired parameters; and matrix A can be computed as

J 11 (, ) () () , . . . , () () TT i ii n i i i n ii      JJ J Aq w q q w J qq w (9)
where is the column vector of the Jacobian matrix for the experiment. is the number of joints.

() i n Jq -th i -th n n
Using usual approach adopted in the identification theory, the estimated unknown parameters can be obtained using the least square method, which yields

 θ ˆ,   X Π k  1 11 ˆ• TT ii i i mm ii            BB X    B p (10)
Using this expression, it can be proved that the covariance matrix for the identification errors in the parameters X can be computed as

1 2 1 ĉov( ) m T ii i         XB B  ( 11
)
where  is the standard deviation (s.t.d.) of the measurement errors. Hence, the impact of the measurement errors on the parameter identification accuracy is defined by the matrix sum that is also called the information matrix.

1 m i   T ii BB
It is obvious that, from practical point of view, the covariance matrix should be as small as possible. However, strict mathematical definition of this notion is not trivial and a number of different approaches are proposed in literature. In most of the related works, the optimal measurement poses are obtained based on minimization of the covariance matrix norm [START_REF] Atkinson | Done Optimal Ex v[END_REF]). This approach may provide a solution, which does not guarantee the best position accuracy for typical manipulator configurations defined by the manufacturing process. Thus, here it is proposed the industry-oriented performance measure, 2 0  , which is defined as the mean square error in the end-effector position after compensation.

To develop this approach, let us introduce several definitions: D1: Plan of experiments is a set of robot configurations and corresponding external loadings W that are used for the measurements of the end-effector displacements and further identification of the desired parameters.

Q

D2. The accuracy of the error compensation 0

 is the distance between the desired end-effector position and its real position achieved after application of error compensation technique.

D3.

The manipulator test-pose is one or set of robot configurations 0 and corresponding external loadings 0 W for which it is required to achieve the best error compensation (i.e.

).

Q mi  2 0 n 
In the frame of the adopted notations, the distance defining the error compensation accuracy can be computed as

0 (     pB X X ) ( 12
)
where the vectors  and are the true parameters values and their estimates, respectively. Matrix 0 B corresponds to the test pose (see expressions (8)). Further, taking into account that

X X  p
 is a function of the unbiased random variables , it can be easily proved that the expectation

 1 ,..., m εε E( ) 0   p
. Besides, the variance can be expressed as

 TT  00 Var( ) E    pX B B X 0 T     (13)
where is the difference between the estimated and true values of the parameters. Expression. ( 13) can be rewritten as 00 t and after relevant transformations in accordance with (10), ( 11), yields the desired expression for the compensation accuracy

   XXX race( E( ) ) TT  BX X B 1 22 00 1 trace ( ) ( ) m T ii i             BB q B qB (14)
As follows from this expression, the proposed performance measure can be treated as the weighted trace of the covariance matrix (11), where the weighting coefficients are computed using the test pose.

Hence, the identification quality (evaluated via the error compensation accuracy) is completely defined by the set of matrices   1 ,..., m BB that depend on the manipulator configurations   1 ,..., m qq . Optimal selection of these configurations will be in the focus of next Subsection.

Problem of the Measurement Poses Selection

Based on the performance measure presented in the previous Subsection, the corresponding optimization problem of the measurement pose selection can be defined as

1 0 , 1 trace ( ) ( ) min subject to ( , ) 0, 1, ii m TT ii i ii i ir             qw BB q B qB Cqw 0   (15)
Here, the matrices describe some constraints, which should be taken into account while solving optimization problem. These constraints are imposed by the work-cell design particularities and usually include the manipulator joint limits, the work-cell space limits, measurement equipment limitations, etc. It should be also mentioned that some constraints are imposed to avoid collisions between the workcell components and the manipulator. Besides, some directions of the applied loading are preferable for the reason of practical implementation.

(, )

ii i

Cqw

It should be mentioned that the component of the matrices i vary with different calibration cases and may include some very specific constraints. For instance, for the case of elastostatic calibration, they can be expressed as

q max i q F is the robot maximum payload,
is the minimum height between the end-point of the calibration tool and the work-cell floor, is the minimum radius to avoid collisions between the applied loading and robot body, min z p min r min  is the minimum angle between the direction of calibration tool and z-axis of robot base frame to ensure that the vertical loading can be applied, and are the boundaries of work-cell space. For the case of geometric calibration, the problem of applying external loading does not exist. So, 2 C , 3 are zero matrices, while , remain the same as in elastostatic calibration.

C

The procedure of solving such an optimization problem could be very tedious for the case when numerous measurement configurations are required for the calibration experiment. For this reason, the problem of interest is to find reasonable number of different measurement configurations and to improve the efficiency of optimization routines employed in the measurement pose selection.

MEASUREMENT POSE SELECTION TECHNIQUES

To solve the above define problem, several techniques can be applied. This section presents the analysis and propose some approaches allowing to obtain acceptable solution in reasonable time. The main difficulties here are related with a large number of variables and complex behaviour of the objective function that has many local minima.

Using Conventional Optimization Techniques

The simplest way to solve this problem is to apply conventional optimization techniques incorporated in commercial mathematical software. It is clear that straightforward search with regular grid is non-acceptable here because of high complexity and enormous number of solutions to be compared. For this reason, three other algorithms have been examined: (i) random search, (ii) gradient search, and (iii) genetic algorithm. Their comparison study is presented below and summarized in Tables 1 and2, where two criteria have been used: computational time and the ability to find optimal solution (evaluated via 0  , the manipulator accuracy after calibration). For all computational experiments, it was assumed that the s.t.d. of the measurement errors is 0.03mm.

The benchmark example deals with the calibration experiments design for 6-dof industrial manipulator KUKA KR270, whose nominal parameters can be found on the manufacturer website (www.kuka.com). The robot has a serial architecture with six actuated revolute joints, so 24 independent geometric parameters should be identified in general case. But for this example, to reduce computational efforts and evaluate the algorithm capability before applying to the problem of real dimension, only nine of the most essential parameters were identified (which have major impact on the positioning accuracy). This allowed us to obtain realistic assessments of the conventional optimization techniques capabilities with respect to the considered problem where the number of design variables is high enough (72 for 12 configurations).

The first of the examined algorithm (i) is based on the straightforward selection of the best solution from the set of ones generated in a random way. For this study, 10,000 solutions were generated for different numbers of measurement configurations . As follows from the obtained results (see Tables 1 and2), this algorithm is very fast and requires less than 2 minutes to find the best solution. However, this solution is essentially worse than the optimal one (by 15-30%).

 3, 4, 6,12 m  

The second algorithm (ii) employs the gradient search with built-in numerical evaluation of the derivatives that is available in Matlab. The starting points were generated randomly and, to avoid convergence to the local minima, the optimization search has been repeated 5000 times (starting from different points). In this case, it has been obtained the best result in terms of the desired objective 0  , but computational cost was very high (it can overcome a hundred of hours). So, this technique is hardly acceptable in practice. It is worth mentioning that reduction of the iteration number is rather dangerous here, because there are a number of local minima that the algorithm can converge to (see Table 1 that includes the minimum, maximum and average values of 0  obtained for random starting points). Moreover, as follows from our experience, 5000 iterations are also not enough here. The third of the examined techniques (iii) applies genetic algorithm (GA) that is based on adaptive heuristic search. The optimization has been carried out for 100 times with population size 50 and 20 generations (initial populations were randomly generated). For illustrative purposes, Fig. 1 presents the efficiency of this algorithm for selection of three optimal measurement configurations. It shows the algorithm convergence as well as divergence of the optimal solutions with respect to computing time. As follows from this figure, the optimization results are highly sensitive to the selection of initial population. In particular, the diversity of the optimal solutions got from sequential GA runs is about 25%. So, to achieve the global minimum, the GA should be repeated many times, which leads to essential increase of the computational efforts (more than 6 hours of computations for the considered example). However, compared to gradient search, GA provides acceptable accuracy (only 2% worse) while the computational time is 4 times less. As follows from the ed re the m obtain sults, rando search is rather fast but inefficient here, since it may produce nonacceptable solutions. In contrast, the gradient search is able to find the global minimum provided that it is repeated many times with different starting points. As a compromise, the GA provides intermediate results in terms of accuracy and computational time. However, for problems of the real industrial size, the performances of the GA are also not sufficient. For this reason, the following Subsections are devoted to the improvement of the numerical optimization techniques employed in selection of optimal measurement poses.

Applying Parallel Computing

Since the considered problem requ of the optimization with different initial values, applying parallel computing looks attractive to speed up the design process and to take advantage of multi-core architecture available in modern computers. To evaluate benefits benchmark example has been considered and two algorithms have been examined: (ii)' parallel gradient search, and (iii)' parallel GA with the same parameter settings. The computations were carried out on the workstation with 12 cores. The obtained results are presented in Table 3, which gives the computational time for different number of measurement poses (the attained value of the objective function 0  is very close to those presented in Table 1).

The obtained results are quite expected and confirm essential algorithms using Number of poses reduction of computational efforts. For both optimization methods, the consumed time has been decreased by the factor of 10-12 (compared to the results in Table 2). However, it is not enough yet to solve the problem of real industrial size, where several dozen of parameters should be identified (instead of nine in the benchmark example). To take the advantages of bot efficiency of the parallel implementation, a hybrid technique has been developed. It should be mentioned that some software packages (Matlab, etc.) already implement this idea and use the final solution from GA as the initial point of gradient search. However, since the randomly generated initial populations in GA may cause high diversity of the optimal solutions, the selection of these initial values is also an important issue. For this reason, the embedded hybrid option in GA cannot be directly used and requires additional modifications.

To improve the efficiency of the existing technique, the starting point selection strategy for the gradient search has been modified. To ensure better convergence to the global minimum, it has been proposed to use the best half of final points obtained from GA as the starting points for the gradient search. From our point of view, it ensures better diversity of the starting points and allows to avoid convergence to the local minima. presents the convergence of the hybrid method for the problem of optimal selection of three measurement configurations studied in the previous Subsection (see Fig. 1). It shows the initial points (obtained from GA), optimal solutions as well as the solution improvement with respect to time. As follows from the figure, the hybrid algorithm can converge much faster, but if the number of measurement poses is increased up to 12, the computational time is over 1.6 hour that is still unacceptable for industry.

New Approach: Reducti

Generally, as follows from the identification theory, th way to improve calibration accuracy is to increase the number of measurements (provided that the reduction of the measurement errors is not possible). However, for the manipulator calibration problem, each measurement is associated with a certain robot configuration that also has influence on the final accuracy. It is clear that the best result is achieved if all measurement poses are different and have been optimized during the calibration experiment planning.

s been tions provides almost the same performance as "fulldimensional" optimal plan. Obviously, this reduction of the measurement pose number is very attractive for the engineering practice.

CONCLUSIONS

On the other hand, as follows from our experiences, the diversity of the measurement poses does not contribute significantly to the accuracy improvement if m is high enough. This allows us to propose an alternative which uses the same measurement configurations several times (allowing to simplify and speed up the measurements). This approach will be further referred as "reduction of problem dimension".

To explain the proposed approach in more details, let us assume that the problem of the optimal pose selection ha solved for the number of configurations that is equal to m , and the obtained calibration plan ensures the positioning ac- It is clear that the calibration accuracy 0 km  for strategy #1 is better than the accuracy corresponding to t e s h trategy #2 that can be expressed as 0 m k  . However, a follows from our study, this difference is not high if m is larger than 3. This allows us to essential e the size of optimization problem employed in the optimal selectio of measurement poses without significant impact on the positioning accuracy.

To demonstrate the validity of the proposed approach, the benchmark example has been solved using strategies # s ly reduc n 1 and #2 assuming that the total number of measurements is equal to 12 (i.e. using different factorizations such as 12 1  , 62  , 43  , 34  ). Relevant results are presented in Table 4 (see the last line). As follows from them, the factor n izatio 12 1  e easurement poses are different is only 6% better compared to the factorization 34  where measuremen repeated 4 times in 3 different configurations. At the same time the factorizations 62  a 3  give almost the same results as the optimal factorizations 12 1  . On the other hand, the computational time of the op l pose generation for 3 m  is much lower than for 1 m wher all m ts are nd 4 tima 2  . This demonstrates the efficiency of the proposed approach and justify its validity. a new technique for o in robot calibration. works, it is proposed to evaluate the quality of calibration plans via the manipulator positioning accuracy for a given test pose, and to reduce number of design variables in the related optimization problem by means of repeating measurements with lower number of configurations. This technique allows us to essentially reduce the computational time for solving the problem of real industrial size. The advantages of the developed technique were confirmed by a simulation example, where the proposed approach permitted to decrease the computing time by more than 10 times while losing only 6% of manipulator positioning accuracy. Future work in this direction will deal with the efficiency improvement of the manipulator elastostatic calibration. ACKNOWLEGMENT like to acknowledge the f ence Nationale de la R ANR-2010-SEGI-003-02-COROUSSO), France.

  all unknown parameters (both geometric and elastostatic ones), and the matrix B varies depending on different calibration cases

Figure 1 .

 1 Figure 1. Efficiency of GA for selection of three measure-

Figure 2 .

 2 Figure 2. Efficiency of the hybrid approach for selection of

  these notations, let us evaluate the calibrtion accuracy for two alternative strategies that employ larger number experiments km : proposed): the measurement poses are obtained by simple repetition the configurations got fro dimensional optimization problem of size m .

   2,

Table 1 .

 1 Efficiency of conventional optimization techniques

	Algorithm		m 	3	Number of poses 4 6 m  m 	m 	12
			min 0	[mm]	0.0825	0.0607	0519 0.	0.0360
	Random Search		mean 0	[mm]	1.9155	0.1939	0.0905	0.0500
			max 0	[m m]	29.878	14.9746	5608 0.	0.0804
			min 0	[mm]	0.0637	0.0521	0426 0.	0.0302
	Gradient Search		mean 0	[mm]	0.0862	0.0610	0.0477	0.0335
			max 0	[m m]	0.6534	0.2333	1319 0.	0.0681
	Genetic Algo-	min 0 mean  0 	[mm] [mm]	0.0638 0.0689	0.0521 0.0529	0427 0. 0.0433	0.0301 0.0305
	rithm		max 0	[mm]	0.0802	0.0615	0446 0.	0.0309

Table 2 .

 2 Computational time of examined algorithms

		Algorithm	m 	3	Number of poses 4 m  6 m 	m 	12
	Random Search	41s	47s	1min	1.7min
	Gradient Search	24.2h	37.5h	56.3h	103.6h
	G	enetic Algorithm	6.5h	8.3h	10.5h	15.4h

Table 3 .

 3 Computational time of examined parallel computing

	Algorithm	m 	3	m 	4	m 	6	m 	12
	Gra	Parallel dient Sear ch	2.1h	3.2h	4.9h	8.9h
	Parallel Genetic Algorithm	36min	41min	52min	1.5h


A Klimchik

Y Wu

G Abba
email: abba@enim.fr

S Garnier

B Furet

A Pashkevich
email: anatol.pashkevich@mines-nantes.fr

Robust algorithm for calibration of robotic manipulator model

Keywords: Robot calibration, robust identification, industrial robot, weighted least square. 

The paper focuses on the robust identification of geometrical and elastostatic parameters of robotic manipulator. The main attention is paid to the efficiency improvement of the identification algorithm. To increase the identification accuracy, it is proposed to apply the weighted least square technique that employs a new algorithm for assigning of the weighting coefficients. The latter allows taking into account variation of the measurement system precision in different directions and throughout the robot workspace. The advantages of the proposed approach are illustrated by an application example that deals with the elasto-static calibration of industrial robot.

INTRODUCTION

The problem of robot calibration is in the focus of the research community for many years (Stone 1987, Hollerbach 1989, Elatta 2004). However, there is a very limited number of works that directly address the issue of the identification accuracy and reduction of the calibration errors (Mooring 1991, Sun 2008, Hollerbach 2008). In general, there exist two main methods to improve the identification accuracy without increasing the number of experiments. The first of them consists in optimization of the manipulator measurement configuration used for the calibration experiments. The second method deals with proper tuning of the identification algorithms used for estimation of the manipulator parameters (fine choosing of the weighting coefficients, for instant). As follows from the literature analysis, the first method has been studied in a number of papers (Khalil 1991, Borm 1991, Daney 2002), while the second one received less attention in robotic research. For this reason, taking into account particularities of the measurement system used in our experiments, this paper focuses on the enhancement of the second method, which looks rather promising here.

To identify the desired parameters, most of the robot calibration procedure employ the ordinary least-square technique, where all identification equations are treated similarly, with the same weights. This approach perfectly suits to the measurement systems that provide roughly the same precision in all directions and in all workspace points. Mathematically, it corresponds to the i.i.d.-hypothesis concerning the measurement noise (i.e. to the assumption that all measurement errors are unbiased, independent and identically distributed). However, in this study, at least one of these assumptions is violated because the precision of the laser-tracker used in the calibration experiments essentially depends on the direction and the target marker location in the manipulator workspace.

To overcome this difficulty, the weighted least-square technique can be applied. As known from literature, for the linear regression it gives rather good improvement and allows essentially reducing the measurement errors impact. It is evident that for the robotic calibration problem, similar benefits can be gained, but the weighting coefficient selection is non-trivial here because of the high non-linearity of the equations describing the robotic manipulator.

To address this problem, the remainder of this paper is organized as follows. Section 2 defines the research problem and also presents some experimental data concerning statistical properties of the measurement errors. In Sections 3, the identification algorithm is presented that is based on the weighted least-square. Section 4 deals with selection of the weighting coefficients that ensure robustness of the identification algorithm with respect to the measurement noise. Section 5 presents an application example, which illustrates benefits of the proposed approach. And finally, Section 6 summarizes the main contributions of the paper.

PROBLEM STATEMENT

In robotics, the calibration procedure can be treated as the best fitting of the experimental data using corresponding manipulator model ,,, ,

    0 arg min ( ,
i ii f   kp L k p Π q Π , ) i F (1)
which describes its geometrical and elastostatic behavior defined by the known function (.) f whose parameters should be tuned. Here i is the vector of measurements (the Cartesian coordinates of the end-effector target points), the p vector L collects all geometrical parameters, i is the vector of actuated coordinates, the vector collects errors in geometrical parameters, the vector collects all compliances of the manipulator, the vectors and i are used for elastostatic calibration only and correspond to the measurement position without loading and applied external loading respectively.

q 0 i p Π k F
In usual engineering practice it is assumed that all measurements ( i and ) are corrupted by the same measurement noise

p 0 i p 2 ) (0,
N  , which induces errors i ε with zero expectation and diagonal covariance matrix ii . However, for many measurement devices such as laser-trackers, the precision highly depends on the measurement direction and vary throughout the robot workspace. In this case, the covariance matrix can be still assumed to be diagonal, but with non-equal principle elements, i.e.

i i xi yi zi () i  ε () T Ed E • 0 2 E  εε I () T  22 (, i a g 2 , )    εε
, where xi

 , yi  , zi 
are different and vary from measurement to measurement. This phenomena is illustrated by experimental data presented in Table 1, which includes dispersions of the measurement errors for the Cartesian coordinates x, y, z for several measurement configurations used in conventional calibration experiments. These results has been obtained by processing the measurement data for 15 configurations and 18 independent experiments for each of them (i.e. by using data set which consists of 810 values). It should be noted that the manipulator repeatability, which is about 60μm, does not have influence on the presented results because of specificity of the measurement experiments, where only difference in the Cartesian coordinate variations were evaluated (before and after applying external loading). As follows from the presented results, the measurement error dispersion vary from 17μm to 153μm and highly depends both on the direction (x, y or z) and the end-effector location in the manipulator workspace (corresponding to the configurations #1-#15). In particular, for the same configuration xi  , yi  , zi  can differ by the factor 5. Besides, from one configuration to another corresponding dispersions may differ by 7 times. Hence, usual assumptions incorporating in robot calibration techniques concerning measurement noise properties should be revised. This motivates the principle goal of the paper that is aimed at developing of the robust identification algorithm that takes into account variations of xi  , yi  , zi  .

IDENTIFICATION ALGORITHM

In manipulator geometric calibration, the basic expression i usually written as follows () (• ,)

gp ii i   Lq pJ Π  (2)
where

g i p () p i J
is the difference between the computed via direct geometrical model and measured end-effector position, the matrix is the geometrical Jacobian and the superscript '(p)' specifies only components that are related to the robot position. It is clearly the linearized model, but it is valid here since in practice the geometrical errors are low enough compared to the nominal values of the manipulator parameters.

Similarly, for the elastostatic calibration, the basic expression can be presented in the form () (,,• )

ep ii i i  Lq F pA k (3)
where

e i
p is the vector of the end-effector displacements under the loading i , the vector collects all compliances of the manipulator, the matrix

F k () p i

A defines the mapping between the joint compliances k and the and-effector displacements. For the further convenience, both expression (2), (3) can be integrated in a single one () ,• (, )

p ii i i  Lq F pB X (4)
where i p is the vector of measurements and

  ,  X Π k
the vector of the unknown parameters that should be identified, () p i B is corresponding matrix function.

Using above defined notations, the calibration can be reduced to the following optimization problem

() () 1 () () m pT p ii ii i F        X BX p BX p m i n (5)
that can be solved using the least square approach. It leads to the following solution

    1 () () () 11 ˆ• T mm pp p ii i i ii     XB B B T  p (6)
If the measurement noise is Gaussian (as it is assumed in conventional calibration techniques) the covariance matrix for the parameter estimates can be computed as follows

X   1 2( ) ( ) 1 ĉov( ) T m pp ii i      XB B T m (7) where () () 1 p p ii i 
 BB is the so-called information matrix.

However, if the measurement noise varies from configuration to configuration the previous expression should be revised. It can be proved that in general case the covariance matrix is expressed as

     11 () () () 2 () () () ĉov( ) • • • • TT T pp p p pp aa a a aa   XBB BΣ BBB  (8)
where the matrix 12 3 ( , ,..., )

m diag    Σ () p i B
describes the statistical properties of the measurement errors and the matrix aggregates all matrices from expression (4).

() p a

B

In order to improve the identification accuracy, it is reasonable to modify the objective function ( 5) and rewrite it in the form

() 2 () 1 () () m pT p iiii i F        X BX p WBX p m i n )  a  p a (9)
where 12 3

is the matrix of weighting coefficients. This leads to slightly different expression for the parameter estimates (, . . , ,.

m w dag w w i  W    1 () 2 () () 2 ˆ•• • •• TT pp p aa a   XB W B B Wp (10)
where a aggregates . In this case, the covariance matrix can be computed as Hence, to assign the weighting coefficients, the measurement noise variance should be known. However, in practice, exact values of xi  , yi  , zi  are unknown and the estimates can be used only. On the other hand, as follows from practical experience, small variations in the weighting coefficients are not critical and they do not affect significantly the identification accuracy. Nevertheless, if the weights are assigned using small number of experiments, the identification results can be unpredictably affected. Therefore, the problem of computing of the weighting coefficients that are able to ensure robustness of the identification algorithm is important here and it will be in the focus of next Section.

p i p     1 () 2 () 1 () 2 2 2 () () 2

ASSIGNING WEIGHTING COEFFICIENTS

It can be proved that for the linear model, the weighing coefficients ensuring the lowest variance of the unknown parameters can be computed as

/ i wa i   ( 12
)
where i  is the standard deviation of the measurement noise for the i-th identification expression and the constant is the scalier factor introduced to avoid the problem of the units. For example, for the linear regression with a single scalar parameter, this approach allows to reduce the variance from

a 22 22 / ( i () ii
Applying this idea to the problem of robot calibration, it is possible to transform the covariance matrix expression (11) to the form

  1 () 2 () ĉov( ) • • T pp aa   XBΣ B (13)
where all elements are essentially lower compared to (8).

Hence the problem of interests is to obtain i  for each experiment (and for each coordinate), which will be used for computing weighting coefficients .

w i      ( 14
)
where  is a scalar factor that allows to tune the impact of i  , 0  is normalization factor that also allows to avoid division by zero.

APPLICATION EXAMPLE

The developed calibration algorithm has been applied to the industrial robot KR-270. To take into account the compensator influence while remaining approach developed for serial robots without compensators (Pashkevich 2011), an equivalent virtual spring with non-linear stiffness depending on the joint variable 2 is used. Using this idea, it is convenient to consider several independent parameters 2i corresponding to each value of 2 . So, the set of desired parameters 21 22 3 6 can be denoted as the vector . This allows us to obtain linear form of the identification equations. q ,...

k k q ( , ...), , kk k k
To find optimal measurement configurations, the design of experiments has been carried out using industry-oriented performance measure proposed in (Klimchik 2012) for five different angles that are distributed between the joint limits. For each 2 three optimal measurement configurations have been found taking into account physical constraints that are related to the joint limits and the possibility to apply the gravity force (work-cell obstacles and safety reasons). The results of the calibration experiment design are presented in Table 2.

2 q q ) x x   to 22 1/ / ii x  

Fig. 2. End-effector used for elastostatic calibration experiments

To generate elastostatic deflections, the gravity forces have been applied to the robot end-effector (see Fig 1) using specific calibration tool (see Fig. 2). To obtain the desired set of initial data, the manipulator sequentially passed through several measurement configurations. Using the laser tracker with the claimed precession 10 m  , the Cartesian coordinates of the reference points have been measured twice, before and after loading. To increase identification accuracy, three reference points (markers) have been used and the loading of the maximum allowed magnitude 250-280 kg have been applied. In addition, to ensure high identification accuracy for each configuration, the experiments were repeated six times. In total the experimental data include 270 measurements which give 810 equations for identification of 9 desired parameters.

The identification has been performed using Ordinary Least-Square (OLS) and Weighted Least Square (WLS) techniques.

For the second approaches the weights have been obtained using non-compensated deflections after identification joint compliances, normalization factor 0  has been assigned to the claimed precession of the laser tracker (10 m  ).

Corresponding values of the elastostatic parameters are presented in Table 3, where for WLS the results for 0.5, 1, 2   are proposed. It also includes the confidence intervals computed as 3  , where the standard deviation  has been evaluated based on the experimental data using expressions ( 8) and ( 13). The results show that the confidence intervals for OLS and WLS have intersections for all parameters of interest, moreover, the confidence intervals for WLS are always inside confidence intervals for OLS and considerably lower (see Table 4). Another important conclusion is that the choice of the coefficient  does not influence on the final results and for simplicity can be assigned to one. It should be noted that the best weighing coefficients can be computed iteratively, where starting from the second iterations the residuals have been computed for the elastostatic parameters identified using WLS with the weights obtained on the previous iteration. This allows us to increase additionally identification accuracy by the factor 3-20 comparing with single iteration WLS. Finally, the identification accuracy have been increased by the factors 27-42 and the errors do not overcome 1.23% for the 6 and is 0.09-0.34%. for the remainder parameters of interest (such non-equivalent distribution of the identification errors is caused by the specific of design of calibration experiment, where the position accuracy after compensation of the elastic deflections has been chosen as a performance measure (Klimchik 2012)). The benefits of WLS are summarized in Table 4. It includes mutual locations of the results for OLS and WLS, identification accuracy and the benefits of new identification algorithm comparing with old one. The convergence of the iterative procedure presented in Figure 3. and Computer-Integrated Manufacturing Journal, vol. 29, issue 2, pp. 385-393. 6. CONCLUSIONS Mooring B., Driels W.M., Roth Z. (1991). Fundamentals of Manipulator Calibration. John Wiley & Sons, NY, 1991. The paper presents a robust identification algorithm for geometrical and elastostatic calibration of robotic manipulator. In contrast to previous works, the proposed algorithm is based on the weighted least-square technique that employs a new strategy for assigning of the weighting coefficients. Such approach allows us to take into account variation of the measurement system precision in different directions and throughout the robot workspace. Application of this technique leads to the essential improvement of the identification accuracy. The advantages of the proposed approach are illustrated by an application example that deals with the elasto-static calibration of an industrial robot which is used for milling of large-dimensional composite parts. It shows that the proposed algorithm allows to reduce most of the identification errors down to 0.3%, which is 30-40 times better compared to the conventional technique. Pashkevich A., Klimchik A., Chablat D. (2011)

  follows from detailed analysis, the best selection of the weighting matrix corresponds to the equation •  W Σ I .

Fig. 1 .

 1 Fig. 1. Experimental setup for the identification of the elastostatic parameters

k

  Figure 3a shows variation in 21 k from iteration to iteration and Figure 3b provides corresponding Confidence intervals. This results shows that parameter value does not change k ACKNOWLEDGMENT significantly (3.3% from the original one), however its precision has been increased by the factor 40. After 10 iterations the values of and CI are stabilized 21 The work presented in this paper was partially funded by the ANR, France (Project ANR-2010-SEGI-003-02-COROUSSO). The authors also thank Fabien Truchet, Guillaume Gallot and Joachim Marais for their great help with the experiments.

  Y., Gen L.P., Zhi F.L., Daoyuan Y., Fei L. (2004). An Overview of Robot Calibration. Information TechnologyJournal, vol. 3, no. 1, pp. 74-78. Hollerbach J. M. (1989). A Survey of Kinematic Calibration.The Robotics Review 1, Cambridge, MIT Press, pp. 207-242. Hollerbach J.M., Khalil W., Gautier M. (2008). Model Identification, Springer handbook of robotics, Part B, pp 321-344 Fig. 3. Convergences of identified value for k 21 and its CI Khalil W., Gautier M., Enguehard Ch. (1991). Identifiable parameters and optimum configurations for robots calibration. Robotica, vol. 9, pp. 63-70. Klimchik A., Wu Y., Pashkevich A., Caro S., Furet BModel Calibration of Anthropomorphic Manipulators. Applied Mechanics and Materials, vol. 162, pp. 161-170 Hence using new identification algorithm the elasto-static parameters of the industrial robot Kuka KR-270 have been identified with a higher accuracy. This allows us to increase the efficiency of the error compensation technique proposed in (Klimchik 2013). Klimchik A., Pashkevich A., Chablat D., Hovland G. (2013). Compliance error compensation technique for parallel robots composed of non-perfect serial chains. Robotics

Table 1 . Dispersions of measurement errors in deflections for different test configurations

 1

	Configuration	σ x , [ m] mean std	σ y , [ m] mean std	σ z , [ m] mean std
	Conf. #1	150 ±1	64 ±1	33 ±1
	Conf. #2	57 ±4	86 ±8	118 ±15
	Conf. #3	97 ±9	70 ±5	44 ±8
	Conf. #4	28 ±1	19 ±1	35 ±1
	Conf. #5	72 ±3	48 ±4	17 ±1
	Conf. #6	153 ±8	46 ±3	22 ±1
	Conf. #7	112 ±6	66 ±3	53 ±4
	Conf. #8	74 ±5	55 ±3	59 ±1
	Conf. #9	80 ±9	63 ±7	102 ±15
	Conf. #10	69 ±2	73 ±1	79 ±1
	Conf. #11	80 ±3	36 ±1	26 ±3
	Conf. #12	53 ±4	39 ±1	29 ±1
	Conf. #13	26 ±1	29 ±1	29 ±1
	Conf. #14	88 ±4	121 ±1	42 ±1
	Conf. #15	90 ±6	52 ±3	50 ±1

Table 2 . Measurement configurations Joint angles, [deg]

 2

	q 1	q 2	q 3	q 4	q 5	q 6
	79.20		-5.57	51.00	-97.52	-91.67
	63.00	-0.01	-12.22	-56.49	41.42	150.55
	63.00		-47.98	-70.04	-61.55	177.16
	95.00		33.00	129.69	-98.10	90.57
	95.00	-25.24	-107.01	109.95	-61.19	174.21
	105.00		14.30	55.21	41.26	-152.97
	56.60		44.54	-55.11	41.90	152.06
	56.60	-56.9	64.73	-129.65	-98.260	-90.55
	144.80		104.49	-69.41	61.67	-6.33
	-41.00		-91.68	55.12	41.53	-152.48
	-143.00	-99.85	-32.64	110.31	-61.47	-6.29
	-143.00		-72.01	129.65	-98.09	90.82
	133.00		147.68	129.64	-97.90	90.99
	-60.00	-140	7.59	-110.09	-61.36	-174.09
	-60.00		-52.00	-124.89	-41.62	27.78
					Tool
		Gravity compensator		F
				Force censor
				Test loading

Table 3 . Identified values of manipulator elasto-static parameters using different approaches, [rad/N× m]

 3

	k i	OLS	=0.5	WLS =1.0	=2.0
	k 21 0.297 ±0.010	0.287 ±0.0003	0.287 ±0.0003	0.287 ±0.0003
	k 22 0.287 ±0.012	0.277 ±0.0004	0.277 ±0.0004	0.277 ±0.0004
	k 23 0.315 ±0.018	0.302 ±0.0005	0.302 ±0.0005	0.302 ±0.0005
	k 24 0.302 ±0.032	0.293 ±0.0010	0.293 ±0.0010	0.293 ±0.0010
	k 25 0.251 ±0.020	0.246 ±0.0007	0.246 ±0.0007	0.246 ±0.0007
	k 3 0.396 ±0.031	0.416 ±0.0011	0.416 ±0.0011	0.416 ±0.0011
	k 4 3.017 ±0.248	2.786 ±0.0071	2.786 ±0.0071	2.786 ±0.0071
	k 5 3.294 ±0.506	3.483 ±0.0120	3.483 ±0.0120	3.483 ±0.0120
	k 6 2.248 ±0.725	2.074 ±0.0267	2.074 ±0.0267	2.074 ±0.0267

Table 4 . Benefits of WLS k i Mutual location of CI CI OLS CI WLS σ OLS /σ WLS

 4 Menq C. (1991). Determination of optimal measurement configurations for robot calibration based on observibility measure. International Journal of RoboticsResearch, vol. 10, no. 1,

	k 21	3.4%	0.09%	40.5
	k 22 k 23	4.2% 5.9%	0.13% 0.18%	33.2 33.9	REFERENCES
	k 24 k 25 k 3	10.7% 7.8% Borm J., 7.8% 0.33% 33.1 0.27% 30.1 0.26% 28.8
	k 4	8.2%	0.25%	35.1
	k 5	15.3%	0.34%	42.0
	k 6	32.2%	1.28%	27.2	Daney D. (2002). Optimal measurement configurations for
					Gough platform calibration. IEEE International
					Conference on Robotics and Automation (ICRA 2002),
					pp. 147-152.

  . Enhanced stiffness modeling of manipulators with passive joints. Mechanism and Machine Theory, vol. 46, issue 5, pp. 662-679 Stone H. W. (1987). Kinematic Modeling, Identification, and Central of Robot Forward Manipulators, Boston, Kluwer. Sun Y., Hollerbach J.M. (2008). Active robot calibration algorithm, IEEE International Conference on Robotics and Automation (ICRA 2008), pp. 1276-1281.