Salim Maakaroun
email: salim.maakaroun@irccyn.ec-nantes.frwisama.khalil@irccyn.ec-nantes.fr

Philippe Chevrel
email: philippe.chevrel@emn.fr

Maxime Gautier
email: maxime.gautier@irccyn.ec-nantes.fr

Wisama Khalil

Modeling and Simulation of a Two wheeled vehicle with suspensions by using Robotic Formalism

Keywords: Modeling, Simulation, Automotive, Robotics, Simulators, Dynamics, Vehicle

Models, simulators and control strategies are required tools for the conception of secure and comfortable vehicles. The aim of this paper is to present a systematic approach to develop models for dynamic vehicle, focusing on a two wheeled vehicles whose body involves six degrees of freedom. The resulting model is sufficiently generic to perform simulation of realistic cornering and accelerating behaviour in various situations. It may be used in the context of motorcycle modeling, but also in various situations (e.g. for control application) as simplified model for 3 or 4 wheeled (tilting) cars. The approach is based on considering the vehicle as a multi-body poly-articulated system and the modeling is carried out using the robotics formalism based on the modified Denavit-Hartenberg geometric description. In that way, the dynamic model is easy to implement and the system can be used for control applications.

INTRODUCTION

Modeling and simulating vehicle dynamics are fundamental tools for vehicles research and development. They allow understanding the dynamics of vehicles and improving the design in order to ensure the major challenge of having safe, comfortable and economic vehicles. Hence, the goal is to build a mathematical model that illustrates significant aspects of the physical dynamics and then facilitate performance analysis and assess design tradeoffs. In the literature, most of the models proposed are developed for control applications (Sharp, 1971). They are centred on motorcycle behaviours, and neglect some essential aspects such as gyroscopic effect on the steering handle bar or pitch motion due to the suspension system (Weir 1978), (Katayama et al, 1985). Lately, some advanced models have been developed using multi body systems. R.S Sharp provided his model by using Autosim software with a description lack on the applied method (Sharp et al, 2001). Later, Cossalter et al developed a model based on Lagrange Formalism that consist on interconnected rigid bodies together with suspensions and other flexible components, supplemented by sophisticated tire and engine models (Cossalter et al, 2002). This formulation uses absolute coordinates that do not depend on the topological structure of the system. However, this technique leads to a complicated model, hard to implement and requires complex numerical algorithm to solve the DAE's (Shabana, 1994). Therefore to model a complex system (Rajamani, 2006), (Kiencke et al, 2000) in 3D motion, we claim that it is preferable to proceed in a systematic geometrical description, based on the modified Denavit Hartenberg parameterization (Khalil et al 1986). This description allows to automatically calculate the symbolic expression of the geometric, kinematic and dynamic models by using a symbolic software package as SYMORO+ (Symbolic Modeling of Robots) (Khalil et al, 1997). This formulation leads to a minimum set of differential equations from where the constraint equations for the mechanical system are automatically eliminated. This paper concentrates on developing a dynamical model for a two wheeled vehicle (called bicycle) by applying recursive methods used in robotics. The approach elaborates systematically the symbolic equations of motion and makes the implementation of the dynamic model easier. This work can be extended for various complex vehicles, such as narrow electric tilting car, specifically Smera Car from Lumeneo (Lumeneo) and (Maakaroun, 2010a, b). The paper is organized as follows: the global method is described in section 2.1 and applied to the bicycle system in section 2.2. A dynamic model is then elaborated using a recursive Newton-Euler based Algorithm (Khalil et al, 1987) in section 3. Finally, Simulations results are illustrated and commented and conclusions are done. The paper ends with a summarize conclusion.

GEOMETRICAL DESCRIPTION OF THE CAR

Robotic representation of a multi body system

The bicycle is considered as a mobile robot which is a treestructured multi body system composed of n bodies (links) where the chassis is the mobile base and the wheels are the terminal links. The links are numbered consecutively from the base to the terminal links. Each body C j is connected to its antecedent C i (i=a(j)) with a joint that represents a translational or rotational degree of freedom and can be elastic or rigid. a(j) denotes the link antecedent to link j, and consequently a(j) < j .A body can be virtual or real; the virtual bodies are introduced to describe joints with multiple degrees of freedom like ball joint or intermediate fixed frames. The frame R i (O i , xi , y i, z i) which is attached to the body Ci is defined as following:

The z i axis is along the axis of joint i, the u j axis is defined as the common normal between z i and z j . The x i axis is along the common normal between z i and one of the succeeding z axis, where link i is the antecedent of link j and the origin O i is the intersection of z i and x i . The homogeneous transformation matrix i T j between two consecutive frames R i and R j is expressed as a function of the following six parameters (Fig. 1): • γ j : angle between x i and u j about z i • b j : distance between x i and u j along z i • α j : angle between z i and z j about u j • d j : distance between z i and z j along u j • θ j : angle between u j and x j about z j • r j : distance between u j and x j along z j

Fig. 1. Geometric parameters

The generalized coordinate of joint j is denoted by q j , it is equal to r j if j is translational and θ j if j is rotational. In (Fig. 1), since x i is taken along u k, the parameters γ k and b k are equal to zero. We define the parameter σ j = 1 if joint j is translational and σ j = 0 if joint j is rotational. If there is no degree of freedom between two frames that are fixed with respect to each other, we take σ j =2. In this case, the time derivative of q j is zero. The description of the bicycle considered as a multi body poly articulated system uses the Denavit and Hartenberg (MDH) notations that are commonly used in robotics (Fig. 4). The chassis motion is described with Euler coordinates while all the other links are described with the generalized Lagrangien coordinates. According to MDH description and SYMORO+, C 0 is the base attached to the ground. The structure is defined as a robot with a mobile base by considering C 1 attached to C 0 via a blocked joint. The inertial parameters of this base are those of C 1 and the speed and the acceleration are then the ones of the chassis described in his own frame. Let R f be a fixed reference frame attached to the ground. The body C 1 with a location ζ (i.e. position & orientation) gives the system posture in the frame R f . The movement of the chassis in this mixed Euler-Lagrangien model is given by:

Application for the model

⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 , , z y x z y x z y x V V V V ω ω ω ω ω ω ω ω
Where V x1 , V y1 and V z1 are respectively the longitudinal, lateral and vertical translational speed of the chassis. Where 1 ω and 1 ω & are the angular velocity and accelerations of the chassis. According to this description, the geometric parameters of the tree structure are shown in table 1 and the bicycle motion is completely described by the vector q of the 11 generalized coordinates:

[ ] T q 1 ξ ξ = ; [ ] 4 11 6 9 3 1 q q q r r = ξ -ξ [1x6]
is the posture of the chassis (position & orientation) -r 3 and r 9 are the length of the suspensions, -q 6 and q 11 are the angular positions of the two wheels with respect to their revolute axis, -q 4 is the steering angle.

Table 1. Geometric Parameters of the modal

DYNAMIC MODEL

Dynamic parameters

For each link there are 14 standard dynamic parameters (Gautier et al, 1990) composed of 10 standard inertial parameters (Table 2): -J j = [XX j XY j XZ j YY j YZ j ZZ j ]: the six coefficients of the inertia matrix of link j given in the frame R j , -MS j = [MX j MY j MZ j ]: the three components of first moment of link j around the origin of the frame j, -M j : the mass of link j

Table 2. Dynamic parameters of the modal

For each actuated joint j, we introduce:

-I aj as the total inertia of the rotor of motor and the drive transmission. -F vj , F sj as the viscous and coulomb friction parameters. For a flexible joint, we define: -K j as the stiffness of the joint j For joint 3 and 9 we add K 3, K 9, F v3 and F v9 to the parameters listed in Table 2. These parameters represent respectively the stiffness of the springs and the dampers of the suspensions.

External Forces

The external forces applied to the bicycle, which have the most significant impact on vehicle dynamics, are the contact forces between the ground and the tires. These external forces can be modeled (Pacejka, 2002), estimated (Canudas, 2003) or measured at the center of the wheels by using dynamometric wheels. Aerodynamic forces also have an effect on the vehicle behavior, particularly at high speed (> 90 Km/h).

Euler-Lagrange Dynamic model

The mixed Euler-Lagrange model is obtained from two recurrences of the algorithm of Newton-Euler in the following way (Khalil, 2002):

The forward recursive equations can be summarized as follows: for j= 1 to n, we calculate the total forces and moments on each link ] [

j j j i j j j j j i i i j i j j j j j i j j j i i i j i j a q a q A a q × + + = + = = ω σ ω ω σ ω ω ω ω
(1) j j j j j j j j j j j j j j j j j j j j j j j j j j j j i j j j j j j

i i i i j i j j j V MS J J M MS V M F a q a q P V A V × + × + = Υ + = × + + Υ + = ) ( ) 2 ( ) ( ω ω ω ω σ (2)
The backward recursive equations can be summarized as follows: for j= n to 1, we calculate the forces and moments exerted on body B j by its antecedent B i .

) ( ) ( ) ( ( ) ( ) ( ) ( ) ( ) ( j s j j s j j s j s j s j s i ej j j j j j j j j i j i j s j s j ej j j j j j f P m A m M m f A f f f F f + + + = = + + = ∑ ∑ (3)
Where:

j ω and j ω & are respectively the angular velocity and the angular acceleration of body j. -F j and M j are respectively the total forces and moments applied on the body j with respect to O j.

j j j j (4) The inverse dynamic model gives the joint torques as a function of the joint coordinates, speeds and accelerations. The joint forces or torques are obtained by projecting f j or m j on the joint axis z j and by taking into account the effects of friction and elasticity as follows:

O i of frame R i in j a(j) σ j γ j b j α j d j θ j r j 1 0 2 0 0 0 0 0 0 2 1 2 0 0 0 L f 0 0 3 2 1 0 0 π 0 π r 3 4 3 0 0 0 0 0 q 4 0 5 4 2 0 0 π/2 0 0 0 6 5 0 0 0 0 0 q 6 π 0 7 5 2 0 0 π/2 0 -R a 8 1 2 0 0 0 -L r 0 0 9 8 1 0 0 π 0 π R 9 10 9 2 0 0 π/2 0 0 0 11 10 0 0 0 0 0 q 1 1 π 0 12 10 2 0 0 π/2 0 -R a j XX XY XZ YY YZ ZZ MX MY MZ M 1 XX 1 XY 1 XZ 1 YY 1 YZ 1 ZZ 1 MX 1 MY 1 MZ 1 M 1 2 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 M 3 4 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 6 XX 6 0 0 YY 6 0 ZZ 6 0 0 0 M 6 7
j j j j e j vj j sj f j e j f j j j T j j a m f + = Γ ) ( σ σ j j j j j off q K q F q sign F + = Γ + = Γ Γ + Γ + & & ) ( (6)
For the body C 1 there is no projection on the joint axis, so the equations of the chassis will be represented by the total forces f 1 and moments m 1 exerted by link 0 on link 1. Thus the NE equations of the chassis are expressed in terms of Euler variables [ ] The inverse Dynamic mod mobile base can be written as: 11x1) is the vector of centrifugal, Cori terms.

1 1 1 1 1 1 , , ω ω & & V . ( 7
( ) ( q H q q A & & + = Γ ) , , , , , , , ( ) ( ) , K F F g f q q q f f q J q v s e t t t e t T t t t t = + -H (
-J (11xp) generalized efforts representing the projection of the external forces and torques on the joint axes.

-p: number of the components of f e -

t t t q q q & & & , ,
are the angular position, ve of joints including the variables of the chassis such as:

[ ] T t q q q r r z y x q ψ ϕ θ = 4 11 6 9 3 [ ] T z y x z y x t q q q r r w w w V V V q 4 11 6 9 3 1 1 1 1 1 1 1 1 1 1 1 1 [ ] T z y x z y x t q q q r r w w w V V V q 4 11 6 9 3 1 1 1 1 1 1 1 1 1 1 1 1
our of the bicycle, by using the dynamic model obtained from the equation of Newton-Euler. The simulator architecture is described in figure 5: ,,,,, (

] ) ( ) , ( [ )] ( [ ) , ,
1 e t T t t t v s e t t t f q J q q H q A K F F g f q q g q - - Γ = Γ = (9)
As we said in Section 3.3, the torques related to the variables of the chassis are zero and the torques of passive joints are zero, we can conclude:

) 1 6 ( 0 ) 6 : 1 ( x = Γ and Γ 0 ) 11 ( =
Moreover, by adding the stiffness and the friction coefficient to the suspensions joint as in section 3.1, the torque of joint 3 and 9 will be equal also to zero. Hence the global vector of torques will be written as follows:

[ ] T 4 11 0 0 0 0 0 0 0 0 0 Γ Γ = Γ
The matrix A can be calculated by the algorithm of Newton-Euler, by noting from the relation (8) that the ith column is equal to Γ:

) 0 , 0 , 0 , 0 , 0 , , 0 , ( ) (:, i t u q f i A = ( u i is t 10)
he unit (11 x 1) vector, whose elements are zero except an be obtained with the the ith element which is equal to 1. The calculation of the vector H c Newton-Euler method, by noting that

H = Γ if: ) , , , , 0 , 0 , , ( K F F g q q f H v s t t & = (11)
The Matrix J can be calculated (12) e is the unit (p x 1) vector, wh the well known and widely used one is Pacejka Magic formula. This model captures in steady-state motion, the tire/road forces and moments, in algebraic equations form with respect to load, longitudinal slip and lateral slip. Among a six components of the tire/road contact wrench, we have considered the principal ones: longitudinal and lateral forces.

-To keep the tires in contact with the ground, we must add two constraints to the dynamical model. Therefore, the vertical velocity and acceleration of the contact tire/road with respect to the reference frame must be equal to zero.

⎤ ⎡ f V (13) 0 ) 3 ( 12 ) 3 ( 7 = Φ = ⎥ ⎥ ⎦ ⎢ ⎢ ⎣ t f q V & t t t f f q dt d q q dt d V V dt d & & & & ) ( ) ( ) 3 ( 12 ) 3 ( 7 Φ + Φ = Φ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ (14)
quation (8) becomes:

Where:

x resents the vector of the efforts transmitted by

And the direct dynamical model for simulation will be:

(15)

First Scenario: longitudinal behaviour

In this scenario, the bicycle is subject to traction torque

E λ T e t T t t t t f q J q q H q q A Φ + + + = Γ ) ( ) , ( ) ( & & & -Φ (2 11) -λ (2x1) rep joints to respect the constraints. ] ) ( ) , ( [ 0 ) ( ) 2 2 ( e t T t t x t t f q J q q H q A q - - Γ ⎥ ⎥ ⎦ ⎢ ⎢ ⎣ Φ Φ = ⎥ ⎦ ⎤ ⎢ ⎣ λ 1 T ⎤ ⎡ -
applied at the rear wheel. The shape of this torque stimulates an accelerated phase then a decelerated one as shown in Fig 6 . The trajectory is a straight line with an initial velocity equal to 10 m/s. When the bicycle accelerates, a transfer load will occur from the front to the rear part of the vehicle. Consequently, the length of the front suspension will increase, and the one of the rear suspension will decrease (Fig 7 .). Hence, a negative pitch angle will appear with respect to its revolute axis which is taken to the left as in

Second Scenario: lateral behaviour

In this scenario, the system is subject to a desired steering torque applied on the steering column to follow the trajectory imposed by a steering desired angle (Fig 10). We consider in this case that the longitudinal velocity is constant and equals to 3 m/s. However, to maintain the stability, the bicycle must tilt into the corner such that the resultant force of the lateral acceleration and the weight of the vehicle is along the vertical axe of the vehicle. The desired tilt angle will be the roll of the bicycle and it will be equal to:

  Our model is composed of 12 bodies (Fig.2& Fig.3) connected by 11 joints: -C 1 is the chassis -C 3 and C 9 are the front and rear suspensions. Their movement is represented by prismatic flexible joints.-C 6 is the rear driving wheel and C 11 is the front steering wheel.-C 4 is the steering column -C 2, C 5, C 7, C 8, C 10, C 12 are virtual bodies fixed to other links by blocked joints.

Fig. 2 .

 2 Fig. 2. Multi body description of the bicycle

Fig. 4 .

 4 Fig. 4. Geometric description of the bicycle

  i j P are respectively the orientation 3x3 matrix and the position vector of the origin

  j) indicates the bodies whose antecedent is body -f ej and m ej are the external forces and moments appl body C j on the environment. When σ j is equal to 2 (a bloc by σ j or j σ are eliminated.

  ) f 1 and m 1 are equal to zero because there is no body el (IDM) of a tree structure with a (8) -A (11x11) is the inertial matrix of the system olis and gravity is the Jacobian matrix and Jf e is the vector

Fig. 5 .

 5 Fig. 5. Simulator architecture

  test scenarios are considered. The following assumptions are used: -The aero dynamical forces are neglec -The behaviour of the ground vehicles depe on the nature of the interaction between the tire and the road. As said in section 3.2, this interaction has been modeled and

Fig. 9 .

 9 Fig. 7. Front and rear suspension elongation

  Fig. 6. Rear Torque

& & & & & & &

& & & & & &

& & & & & &

& &

& & & & & =

& & & & & & & & & & & & & & & & & =

& & & &

Preprint submitted to 18th IFAC World Congress. Received October 12, 2010.

& & &

ational Journal of Vehicle Mechanics and Mobility, 17(4), 211-229. Khalil, W., Kleinfinger, J. (1986). "A new geometric notation for open and closed-loop robots," Robotics and Automation. Proceedings. 1986IEEE International Confe rence on , vol.3, no., pp. 1174-1179. Khalil, W. and Kleinfinger, J.(1987) "Minimum operations and minimum parameters of the dynamic model of tree Fig. 10. Desired steering angle structure robots,"IEEE J. Robot.Autom.,no. 6,. Khalil, W. and Creusot, D. (1997) (1971). "The stability and control of motor-6. CONCLUSIONS cycles," Journal Mechanical Engineering Science, Based on a mixed Euler-Lagrange formulation, a two wheeled vehicle dynamic model was developed for simulation's purpose and some dimensioning and control applications. Future works will consist on applying this technique on more complex and sophisticated vehicles as narrow tilting cars such as the Smera from our collaborator Lumeneo. Two approaches will be compared : i) a bottom up one consisting to analyze the ability of the simplified model proposed in the present paper to give realistic results with ad hoc parameters using identification from experimental data, ii) a top down one, showing how to simplify the complete model for control purpose.. vol. 13, no. 5, pp. 316 -326. Sharp, R. S., (2001) x, y, z: longitudinal, lateral and vertical distance of the centre of gravity with respect to the reference frame. θ, φ, ψ: roll, pitch and yaw of the bicycle Ra: Radius of the wheels: 0.275 m REFERENCES h: 0.625 m :height of c.g of bicycle from ground Lf = Lr = 0.85 m: longitudinal distance from c.g to front and rear wheels; g = 9.81 m/s 2 : gravitational constant Canudas deWit, C., Tsiotras, P. and al, (2003)."Dynamic friction models for road/tire longitudinal interaction," M 1 = 1000 Kg: mass of the chassis Veh. Syst. Dyn.,vol. 39,no. 3, M 3 = M 9 = 1.32 Kg: mass of a suspension Cossalter, V., Lot, R. (2002). "A motorcycle multi-body M 6 = M 11 = 10 Kg: mass of a wheel model for real time simulations based on the natural XX 1 =43.3 Kg m 2 : roll moment of inertia of bicycle coordinates approach," Vehicle System Dynamics, vol.

YY 1 =254 Kg m 2 : pitch moment of inertia of bicycle 37, no. 6, pp. 423-447.

ZZ 1 =270 Kg m 2 : yaw moment of inertia of bicycle Gautier, M. and Khalil, W. (1990), "direct calculations of XY 1 = XZ 1 = YZ 1 = MX 1 = MY 1 = MZ 1 =0 minimum set of inertial parameters of serial robots," XX 6 = YY 6 = XX 11 = YY 11 = 0.207 Kg m 2 IEEE Trans. On Robotics and Automation, Vol. RA-ZZ 6 = ZZ 11 = 0.378 Kg m 2 6(3), p. 368-373.

K 3 = K 9 =50000 N/m: stiffness of a suspension Rajamani R. (2006). vehicle dynamics and control, Springer.

F v3 = K v9 =5200 Ns/m: damping coefficient of a suspension Katayama, T., Aoki, A. (1988). "Control Behaviour of


Paul Chiambaretto
email: paul.chiambaretto@polytechnique.edu

André De Palma
email: andre.depalma@ens-cachan.fr

Stef Proost
email: stef.proost@econ.kuleuven.be

The impact of interregional and intraregional transportation costs on industrial location and efficient transport policies

Keywords: JEL Classification : R11, R12, R42, R48, R58 Economic geography, Industrial location, Transportation costs, Intraregional, Interregional, Concentration

Almost all models of the (New) Economic Geography have focused on interregional transportation costs to understand industrial location, considering regions as dots without intraregional transportation costs. We introduce a distinction between interregional and intraregional transportation costs. This allows assessing more precisely the effects of different types of transport policies. Focusing on two regions (a core and a periphery), we show that improving the quality of the interregional infrastructure, or of the intraregional infrastructure in the core region, leads to an increased concentration of activity in the core region. However, if we reduce intraregional transportation costs in the periphery, some firms transfer from the core to the periphery. From an efficiency point of view, we observe that, in absence of regulation, the concentration of firms is too high in the center. We show what set of policies improves the equilibrium.

Introduction

At the end of the nineteenth century, Marshall (1890) explained that "a lowering of tariffs, or of freights for the transport of goods, tends to make each locality buy more largely from a distance what it requires; and thus tends to concentrate particular industries in special localities." Indeed, during the industrial revolution, Marshall witnessed a key moment in the history of geography. As [START_REF] Bairoch | Regional economics: a new economic geography perspective[END_REF][START_REF] Bairoch | Regional economics: a new economic geography perspective[END_REF] explains, throughout the nineteenth century, transportation costs have been divided by ten, and at the very same time, inequalities between countries emerged: the standard deviation of GDP per capita in Europe has been multiplied by 7.5. A very detailed analysis of this phenomenon is given by [START_REF] Lafourcade | New economic geography : The role of transport costs[END_REF].

The reduction of trade costs is one of the causes of these inequalities.

In this paper, we distinguish between interregional and intraregional transportation costs. This allows assessing more precisely the effects of different transport policies on the regional distribution of activities.

Even if Marshall had a good intuition of the relation between reduction of transportation costs and concentration of activities, the main contributions began to emerge in the 80's. Three main models have been developed to study the effect of interregional trade on industrial location. The first one has been developed by [START_REF] Helpman | Scale economies, product differentiation and the pattern of trade[END_REF] (but it is usually called the Dixit-Stiglitz-Krugman model). He sets the basis of interregional trade in an imperfect competition framework with increasing returns to scale. He considers two sectors (agriculture and industry) that are both present in two regions. A manufactured good can be sold in the two regions, but to sell it "abroad", a transportation cost must be paid by the consumer. The author concludes that a decrease of the transport costs will affect positively the two regions (by diminishing the price level in both regions). Nevertheless, the impact will be higher for the big region than for the small one, and so there will be an increase in inequalities between regions.

This model has been further developed by [START_REF] Helpman | Scale economies, product differentiation and the pattern of trade[END_REF]. Consider 2 regions A and B, and 2 production factors: labor is homogenous and can be used either in the agricultural or the industrial sectors. Labor is immobile, whereas capital is mobile between the two regions. Each worker has a unit of capital he can invest either in region A or B. But the returns he gets from the capital are spent in his own region. Defining equilibrium as the equality of returns in the two regions, they managed to highlight what they referred to as the "Home-Market Effect". This means that the share of industry in the employment of the central region is bigger than its share in the population. The reason is that, in presence of increasing returns and transportation costs, firms will prefer to locate near the biggest market. This effect will be even more intense when trade costs are reduced.

The last main contribution is by [START_REF] Krugman | Increasing returns and economic geography[END_REF]. He analyzes which workers move between regions, and how wages are set endogenously. Workers and firms move between regions comparing their expected utilities and profits. Krugman observes that, below a certain threshold, the reduction of transportation costs will automatically lead to the concentration of all the industrial activity in one of the two regions.

These three models have two points in common: first, they all find that the reduction of interregional transportation costs will increase inequalities between regions. Second, they all use the same assumption: they consider regions as dots, without dimension, in which there are no intraregional transportation costs (see Behrens and Thisse, 2007). This last point is disturbing since a whole field of economics has focused on cities: urban economics. Many contributions have been made, but most of them were neglected by interregional trade economists. The first attempt to unify this field was the paper of [START_REF] Tabuchi | Urban Agglomeration and Dispersion: A synthesis of Alonso and Krugman[END_REF], in which the author proposes a synthesis of [START_REF] Alonso | Location and Land Use[END_REF] and [START_REF] Krugman | Increasing returns and economic geography[END_REF].

Other papers have contributed to the linkage of these two growing fields. We can think of the paper by [START_REF] Duckworth | The Rise and Fall of Regional Inequalities[END_REF],where he observes that, with congestion costs, the "tomahawk curve" of [START_REF] Krugman | Increasing returns and economic geography[END_REF] becomes a bell-shaped curve The idea is that the congestion costs will reduce the incentive for firms to remain in the center when interregional trade costs are reduced. Cavailhès and al. (2006) go one step further and investigate the way the structure of cities is affected by interregional trade, shifting from a monocentric to a polycentric configuration to reduce congestion costs. However, all these contributions have neglected a very important scale: "the region". As observe Behrens and Thisse (2007), these contributions have gone from the interregional scale to the urban one, skipping the region. We help filling this gap, making the difference between interregional and intraregional infrastructures. Such a distinction has already been made by [START_REF] Macmillan | Industrial location and public infrastructure[END_REF]. They use this distinction to analyze FDIs (Foreign Direct Investments) in developing countries. They observe that an improvement of the international infrastructure will motivate firms to move to developed countries, whereas an improvement of the regional infrastructure in the periphery will lead to a transfer of firms from the developed country to the developing one. We want to use such a distinction, between interregional and intraregional transportation costs, to assess more precisely the effects of different transport policies.

There are several examples of the effects of infrastructure on regional inequalities. [START_REF] Vickerman | Infrastructure and regional development[END_REF][START_REF] Pion | Accessibility and Economic development in Europe[END_REF] has studied the effects of the European transport policy on regional inequalities.

An historical example is provided by [START_REF] Cohen | La mondialisation et ses ennemis[END_REF]. He explains that during the French colonization of Algeria, many roads were built to connect distant villages to central cities.

These roads allowed firms from the center to sell their products to the villages. Far from improving the situation, these roads emptied the remote villages and increased the spatial polarization of activities.

It is crucial to understand the various mechanisms at stake, in order to implement a transport policy taking into account the cohesion objective. Another interesting objective is the normative study of such a situation. Indeed, is the geographical equilibrium an optimum from a Pareto point of view? Few contributions have focused on the normative point of view of the economic geography, and none of them made a distinction between interregional and intraregional transportation costs. This paper has two sections. In the first section, we start from an existing model of intraregional trade (the Helpman and Krugman's one), and introduce intraregional transportation costs. With these intraregional costs, we can assess the effects of different transport policies on industrial location. In a second section, we take a normative point of view. We compare the spatial equilibrium to the Pareto optimal one. We show that the geographical equilibrium is not optimal. We next examine policies like road pricing that improve the efficiency of the equilibrium. A numerical example with transport policies for African countries illustrates the model. In a final part of the paper, we draw some conclusions from our model. Detailed mathematical proofs of the propositions are relegated in the appendixes.

A theoretical model to study the impact of interregional and intraregional transportation costs on industrial location

Description of the model

Most previous models in Economic Geography neglected the intraregional transportation costs. One of the exceptions is the model developed by [START_REF] Macmillan | Industrial location and public infrastructure[END_REF], in which they make a distinction between international and domestic trade costs. If their paper is a first step in the understanding of the impact of the intraregional trade costs, they do not use a normative framework: they do not look for efficiency improving policies.

In this paper, we add to the Helpman and Krugman's model (1985), intraregional transportation costs. The advantage of this model is that it has an analytical solution, what it is not the case of the [START_REF] Krugman | Increasing returns and economic geography[END_REF] model. Our reasoning is very close to the one proposed by [START_REF] Macmillan | Industrial location and public infrastructure[END_REF]. We change the mathematical notations so as to obtain results more directly.

We consider 2 regions: A (the center) and B (the periphery), where A has a bigger share of the population. In each region, there are two sub-regions: factories on one side and houses (with shops) on the other side.

Figure 1. Structure of the regions and transportation costs

The structure of these regions allows us to have intraregional as well as interregional transportation costs.

Main Assumptions

Before solving the equilibrium, we need some assumptions. We consider a model of two regions: the center and the periphery. In the economy, there are L workers. A share  of the workers are in the central region, with 1/ 2  

. In this model there are two kinds of factors:

labor which is immobile and capital which is mobile, all members of the population own one unit of capital. As in [START_REF] Helpman | Scale economies, product differentiation and the pattern of trade[END_REF], wages in the two regions are set to 1:

1 AB ww  .
More precisely, we define the different unit transportation costs within regions (AA, BB) and between regions (AB, BA). The intraregional transportation costs are given by AA  and BB  , whereas interregional ones are given by AB BA   .

We suppose that the intraregional transportation costs are more important in the poor region (given the low quality of transport infrastructure) than in the rich region. Moreover, we assume that interregional transportation costs are much higher than intraregional ones. We then have the following inequality:

.

BA AB BB AA      

Equilibrium

The objective of this model is to determine the value of  , which is the share of the industry located in the central region. As labor is immobile, intraregional commuting costs do not affect the location of production. We focus then on transportation costs. To obtain the value  we must first specify some elements concerning production and demand.

Production

The cost function of a firm is defined as    

C q f r cq   , so that there are increasing returns. Each firm needs f units of capital and each of the L workers have a unit of capital.

All firms have the same cost function and produce each a different variety. Denote by  the share of the capital invested in A, Since capital is perfectly mobile between the two regions, the number of firms in regions A and B are:

(1)

A L n f   and   1 B L n f    .
Given we have assumed that the capital is mobile, the spatial equilibrium is defined by the equalization of returns :

      * * * AB rrr    .
The returns are spent in the region of the owner. We obtain the value of the income of each region:

(2)

      1 1 1 AB Y r L and Y r L                  .
Without iceberg transportation cost, the profit equation for a representative firm i is given by

    () i i i i p q p fr cq p    
. Profit maximization leads to the equilibrium price

1 1 i i pc     
, where

ii i ii qp pq     .
With the utility function defined in the Section 2.3.2, we have:

i   , where  is the elasticity of substitution between varieties. Now, we introduce iceberg transportation costs as the percentage of the good that is lost because of transportation costs. If we want to receive a quantity q of a product, it will be necessary to ship q  of the product, with 1   . We then infer the price lm p , paid by a consumer living in the region m, and purchasing a product made in region l :

(3)

1 lm lm lm i c pp      where , l A B  and , m A B  ,
where lm  represents the iceberg cost.

Demand

As in the Dixit-Stiglitz-Krugman model, we have the following utility function:

1 U M A   
where  , the part of the income spent on the composite good, is such that 1   . M is the composite good produced by the industry and A is the numeraire good. This numeraire good can be transported without cost between regions and every unit of labor that is not used to produce the industrial goods can generate one unit of the numeraire good. This allows us to set the price of the numeraire good and the wage equal to 1 in both regions, as in Krugman's model. The share of revenue dedicated to the composite good is

m Y  , where m Y is the total income in region m. Maximizing ( ( ) ) m lm lm i M Y q i p       L
, we obtain the demand for a variety i, made in l and consumed in m :

(4)

      1 ( ) lm i lm m lm j j p q i Y p          ,where , l A B  and , m A B  .
In the rest of this model, we will use the notation:

( 1)

   σ lm lm τ .
We observe that lm  takes a value between 0 and 1, if 1   . When lm  is near 1, there are low barriers for trade.

We add some assumptions on the values of  :

(5)

AA BB AB BA      and (6) 1 AB BB     .
The first hypothesis is just the transportation cost inequality written in terms of  . The second assumption will be justified soon, we need it to make sure the share of industry in the region belongs to [0,1]. In other words, this inequality means that interregional transportation costs must be much higher than the intraregional ones. This assumption seems reasonable.

The total demand for the variety i produced in A is given by the sum of the demand for this variety by the region A and by the region B. The revenues in the equation ( 4) refer to the equation (2). Since the prices are given by the equation ( 3) and the number of firms by (1), we find:

(7)                  1 1 1 1 . 11 AA A AB B i A AA AB AB BB r L r L q c L L L L                              
The first part is the demand for this variety from the consumers of the region A, whereas the second part is the demand from the consumers of region B

Determining the equilibrium

In the long term, the profits are just high enough to cover the cost of the capital. So the profit of the firm, producing the variety i, in the region A is:

(8)                 * * Π 0 A AA AA AB AB AA AA AB AB i p i q i p i q i c q i q i fr            Let define the aggregate production as       A AA AA AB AB q i q i q i  .
Then we have:

(9)   () , ( 1) A A cq i r f    
which can, using the demand functions defined by the equation ( 7), also be written as:

(10)                1 1 1 . 11 AA A AB B A AA AB AB BB rr r f                              Symmetrically, we find () B r  : (11)                1 1 1 . 11 AB A BB B B AA AB AB BB rr r f                             
The spatial equilibrium is obtained when the returns in the two zones are the same. Therefore we are looking for the  value such that

            1 1 . 1 1 1 1 AB BB AA AB AA AB AB BB AA AB AB BB                                     
After simplifications, we find  (the share of industry in the central region):

(12)    ,          AA AB AB BB where      Ψ1 BB AB AB AB AA BB              .
Proposition 1: If the interregional transportation cost is sufficiently large, compared to the intraregional cost (

  1     BB AB )
, there exists a unique interior equilibrium where the industrial activity is shared between the two regions (

[0,1]  ). If the interregional cost is too low (   1     BB AB )
, there is a corner solution and all the industrial activity is in the center ( 1

  ).
The proof of this proposition is relegated to Appendix 1.

2.4.Comparative statics

One of the advantages of building our model on [START_REF] Helpman | Scale economies, product differentiation and the pattern of trade[END_REF] is that we have an analytical solution for the equilibrium. This allows comparative statics. Indeed, we want to know the effects of improving the different types of infrastructure on the industrial location patterns.

In this part, we modify the quality of a specific type of road and we evaluate its impact on the distribution of industrial activity. We assume that the funds (and resources) for the realization of this infrastructure are external. For instance, the funds could come from an international development agency (World Bank) or be part of a federal effort to help peripheral regions (Regional investment Fund in the EU). We show in Appendix 2 the following proposition.

Proposition 2: Both the improvement of the quality of interregional infrastructure and the reduction of intraregional transportation costs in the center lead to a higher concentration of industries in the center. However, lowering the peripheral intraregional transportation costs increases the attractiveness of the periphery for firms.

Proposition 2 can be explained using the notion of market potential. As it was defined by [START_REF] Grasset | The market as a factor in the localization of industry in the United States[END_REF], and then extended by [START_REF] Head | Market potential and the location of Japanese firms in the European Union[END_REF], the market potential is like a weighted sum of the different potential sales on the national market and the surrounding markets where a firm would like to export its product. The weights are inversely proportional to the trade costs.

Since the center is bigger than the periphery, it has naturally a higher market potential, which attracts firms. Reducing the interregional transportation costs will increase the market potential of the center, because central firms will have a better access to the periphery.

Because of increasing returns, firms wish to move to the center. The reasoning is the same with intraregional transportation cost in the center. However, a decrease of the intraregional transportation cost in the periphery will increase the weight of the periphery in its market potential. The market potential of the periphery will then increase, and will consequently attract firms that will relocate from the center to the periphery.

Is this equilibrium efficient and can we improve it?

In this section, we first look for the efficiency benchmark: what geographical distribution of industrial activity maximizes the sum of utilities in the two regions. The equilibrium we described in the previous section is then compared to the benchmark. Next we look for policies that could bring the equilibrium closer to the efficiency benchmark.

Computation of the indirect utility functions in the two regions

In order to maximize overall efficiency, we need an analytical expression for the utility in both regions. Recall that

1 U M A   
where 1   . Moreover, since A is the numeraire, then the budget constraint reduces to: .1 P M A  . Substituting into the utility function, we obtain:

1 1 . () U M PM  



Maximizing the utility, with respect to M, and injecting the optimal value of M in the utility function, we obtain the indirect utility function:

(13)     1 1 . V P       

Moreover, we know that maximizing

V is equivalent to maximizing any increasing transformation of V . We will then rely on the following indirect utilities in this section: (14) and .

A A B B V P V P   
In the new economic geography literature, the price index in the region A is given by:

(15) 1 1 ( . ) 1 A AA A AB B c P n n        
Using equation ( 1) for the number of firms, we obtain equation ( 16), which can be rewritten as equation ( 17).

(16)

1 1

(1 ) , 1

A AA AB c L L P ff               or (17)     , 1 A A A AB K P        where 1 1 1 cL K f           and 1 . 1     
We then have the indirect utility functions for an agent living in A (

A V ) and in B ( B V ) : (18) ( ( 1

))

.

A AA AB VK           
Similarly, we obtain:

(19) ( (1 )) . B AB BB VK           

Maximization of the total welfare in the economy

We start by looking for the maximum of the unweighted sum of utilities in the two regions.

There are two justifications for this. First one can allow for lump sum redistributive transfers by a federal government. Second one can see the maximization of utilities as the basis for an efficient bargaining between the two regions where the two regions share the gains of a better equilibrium via transfers among themselves.

Since a share  of the population is in the center (and 1   in the periphery), the welfare function in the economy is then given by:

  1 , AB W V V               1 1 1 . AA AB AB BB W K K                                  
We look for the value o  that maximizes the welfare function:

            max max 1 1 1 . AA AB AB BB W W                                  
The first order condition gives us the following value:

(20) , o BB AB BB AB AA AB             where      1 1 . 1 BB AB AA AB                 
We can note that 11     , and as a consequence:

[0,1] o  
. The proof of these properties is relegated to Appendix 3.

Comparison of the optimal and equilibrium values of the shares in industrial activity

The question here is whether the concentration in the center is too high when the regulator doesn't intervene. To answer this question, we must compare the equilibrium value Eq  and the optimal value o  . Recall that:

,

Eq BB AB AB BB AB AA AB           
and that :

    . o BB AB AA AB BB AB AA AB BB AB                     
After some calculations (see Appendix 5), it can be shown that Eq o   . The results are summarized in:

Proposition 3: At equilibrium, the spatial concentration of industrial activity in the center is too high compared to the first best optimum.

That the spatial equilibrium leads to a higher concentration than the first best optimum can easily be understood. In fact, without any intervention, there is a kind of "snowball effect"

(the expression was first used by [START_REF]Economic Theory and Underdeveloped Regions[END_REF]. For a given level of infrastructures, the center will attract firms. With these new firms, the market potential of the central region will increase, attracting more firms and so on… If nothing is done by the regulator, there will be a concentration of firms in the center that will be higher than the first-best optimum.

What would become of Prop. 3 when one gives a higher weight to peripheral citizens and there would be no lump sum redistribution possible or no efficient bargaining between the two regions? In Appendix 4 it is shown that a higher weight to peripheral citizens decreases the value of the optimal share of industry in the center.. This is obvious since, giving a higher importance to peripheral citizens, the concentration of firms in the center reduces the price index in the periphery.

In order to improve total welfare, the concentration of firms in the center should be reduced.

This means that the government must intervene so as to favor the transfer of firms from the center to the periphery.

Policies to reach the optimal location : the role of road pricing

In order to decentralize the social optimum o  , we will use a set of incentives. In principle one could use different instruments; the regulator could tax the firms in the center and/or subsidize the firms in the periphery. Here, we focus on instruments that tax or subsidize the use of transport infrastructures. This can be understood as a form of road pricing or as a shadow cost used to compute the optimal size of different transport infrastructures. These incentives will allow us to match o  and Eq  .

The intuition is clear from the outset, we will have to tax the use of the interregional road, tax the use of the intraregional road in the center and/or subsidy the use of the intraregional road in the periphery. The question is then to know the precise value of these taxes or subsidies.

Taxation of the use of the interregional road

We want to tax the interregional infrastructure. This is equivalent to reducing the value of , AB  which is the "freeness of trade". To do this, we will look for the value t such that AB AB AB t   . We want to reach the value o  , so we look for t , that solves :

        ) . 1 ( AB AB AB BB AA BB o AB AB AA BB                   
After some calculations, and solving the polynomial function in AB  , we find :

(21)   2 2 1 ( 1) (1 ) 4 ( )( 1) 2( 1) o o o o o o AB AA BB BB AA AA AB o                                 

Two remarks are in order:

Remark 1: From AB  , we can find the value of .

Taxation of the use of intraregional infrastructure in the center

We want to tax the use of intraregional infrastructure in the center. This is equivalent to reducing the value of AA  , which is the "freeness of trade". To do this, we will look for the value t such that AA AA AA t   . We want to reach the value o  , so we look for t , that solves:

         1 . AA BB AB AB AB BB o AA AB AB BB                   
After some calculations, and solving the polynomial function in AA  , we find (22) (1 )( ) . ( )

BB AB AB AA AB o AB B B B B             
And the above remarks become:

Remark 1: From AA  , we can deduct the value of .

Subsiding of the use of the intraregional road in the periphery

We want to subsidize the use of the intraregional infrastructure in the periphery. This is equivalent to increasing the value of BB  , which is the "freeness of trade". To do this, we will look for the value of the subsidy s such that .

BB B B B B s   
We want to reach the value o  , so we look for ' BB  , that solves:

        1 . () BB BB AB AB AB AA o BB AA AB AB                   
After some calculations, we obtain: ( 23)

2 ( 1) . ( 1) ( )

oo AA AB AB oo AB AA BB                    Remark 1 : From BB  , we can deduct the value of BB s . Remark 2 : Increasing the value of BB  to BB s   is equivalent to reducing the iceberg cost from BB  to BB k   .
The value of k is then given by the following expression:

  1 1 . BB BB BB ks      
We have developed a set of incentives such that the equilibrium will be the optimal one. To do this, we focused on road pricing. One can reach the optimum in three different ways:

taxing the use of the interregional infrastructure; taxing the use of the intraregional infrastructure in the center or subsidizing the use of the intraregional infrastructure in the periphery. Note that subsidizing the intraregional infrastructure in the periphery does not mean building a road. Indeed, the creation of the infrastructure is costly, while a subsidy is in principle a mere transfer of resources that corrects incentives and is not consuming real resources (except for the transaction costs). We illustrate our model below with a numerical example.

A numerical example: Mozambique and Malawi

In order to illustrate the mechanisms of the model, we take a numerical example. Since comprehensive data on interregional and intraregional transportation costs for road transport are difficult to gather, we use from UNCTAD (2004). In their report, they calculate very accurate values for transportation costs in Africa. We take the example of two countries, instead of two regions. This choice is mainly due to the availability of data, but especially because it does not change anything concerning the assumptions and results. The two selected countries are Mozambique (the center) and Malawi (the periphery).

First, we must calibrate our model. To do so, we match real data with the various parameters we use in the model. Table 1a displays values for the central country, whereas table 1b does it for the periphery. Note that several values have been tested for the elasticity of substitution, and they do not change dramatically the results. We can now compute the spatial equilibrium and the spatial optimum. It is interesting to note that the simulated spatial equilibrium ( 0.75

Table 1a. Numerical values for the parameters of Mozambique

 

) is very close to the real value ( 0.74

 

). As predicted by Proposition 3, there is a significant difference between the spatial equilibrium and the optimal concentration. The next step is to calculate the values of the different taxes or subsidies to reach the optimal concentration. From these new values of transportation costs, we can deduct the impact on the product prices: Intraregional in the periphery 4% 2%

This numerical example illustrates the different values predicted by our model. The values of the various taxes or subsidies are not unrealistic. These numerical simulations confirm that the predictions of our model are in line with the theoretical results. At equilibrium, the concentration of activities in the center is too high and the use of interregional road taxes can restore the optimum.

CONCLUSION

If a whole literature has studied the impact of transportation costs on the geographical distribution of activities, this paper is among the first to analyze the specific roles of interregional and intraregional transport costs .

Several conclusions can be drawn from this paper. First, the improvements of the different categories (interregional/intraregional) of infrastructures have different effects on industrial location. We confirm Martin and Rogers (1995analysis: if the decision maker whishes to use the transport policy to reduce regional inequalities, then s/he must improve the quality of the peripheral intraregional infrastructure. The second conclusion follows from the normative analysis. Indeed, we show that the spatial equilibrium is far from being Pareto optimal.

Without any intervention, the spatial equilibrium will lead to a concentration of firms in the center that is too high. The third result is that it is possible to use tax and subsidy instruments such as road pricing to reach the optimal share of firms in the center.

Yet, this model could be improved in several ways. First, it was designed in such a way that interregional and intraregional transportation costs were independent. In reality, mattters are more complex and intraregional infrastructures may affect directly the interregional transportation costs. These network effects must be taken into account if we want a more realistic view of the effects of infrastructures on industrial location. Second, the introduction of congestion costs within regions would probably give interesting results. As it was explained by [START_REF] Lafourcade | New economic geography : The role of transport costs[END_REF], congestion costs moderate the spatial polarization of activities. In our model, they would affect strongly the optimal taxes and subsidies. With these improvements, future research allows to understand more properly the dynamics of spatial activities.

Appendix 1. Proof of Proposition 1

We want to prove that

  1 AB BB     , then [0,1]   Proof of 0   We know that:    0 AA AB AB BB        . If 0
  , then we must have Ψ0  . Let's prove it by contradiction. If Ψ0  , then we have

     Ψ 1 0 BB AB AB AB AA BB                    . 1 BB AB AB AA AB BB              Since 1   and BB AB AA AB        and since AB BB
  , the previous line ca not be true.

This means that Ψ0  , and then 0   .

Proof of 1  

We know that    0

AA AB AB BB        . If 1   , then we must have :         1 BB AB AB AB AA BB AA AB AB BB                          1 0. BB AB AB AB AA BB AB BB                  Since    1 0 BB AB       and   0 AB AA   , we must have 0 BB AB BB       which is true if   1 . . AB BB     Proof that if   1 AB BB     , then 1   We know that:    0 AA AB AB BB        . If 1   , then we must have:         1 BB AB AB AB AA BB AA AB AB BB                         . BB AB AA AB AA AB BB              Since   1 AB BB     , then   1 AB BB       , so we have    ( ). BB AB AA AB BB AA AB           
Moreover, we can note that:

  AA AB AA AB        
. So we can conclude that we have:

     BB AB AA AB AA AB BB             , which means that if   1 AB BB     , then 1   .

Appendix 2. Proof of Proposition 2

Impact of  AB : Improving the quality of the interregional infrastructure is like an increase in

AB  .          2 1 2 Ψ[ 2 ] BB AB BB AA AB AB BB AA BB AB AB AA AB AB BB                                  We have: 0     AB λ if the condition 2 AA AB    
is respected (which always holds with our assumptions). Using this hypothesis, we observe that the improvement of the infrastructure between regions will strengthen the concentration in the center.

Impact of 

AA : We measure the effects of an improvement of infrastructures in the center. We expect that it will lead to a higher concentration in the center.

      2 [] Ψ BB AA AB AB BB AB BB AA AA AB AB BB                       We obtain: 0     AA λ
. As predicted, we conclude that the higher quality of infrastructure in the center will increase the concentration.

Impact of 

BB : Since the other actions on infrastructures have led to a higher concentration, we expect that the reduction of transport costs in the periphery will lead to a reduction of the concentration in the center.

 

        2 1 [ ]Ψ AB AB AA AA AB AB BB AA AB BB AA AB AB BB                                 We find that 0     BB λ if we have AB AA
   (which holds with our hypotheses) .The better quality of infrastructure in the periphery will lead to a relocation of firms from the center to the periphery.

Appendix 3. Remarks concerning the value of ζ

First, we want to show that 01  

We can rewrite  as:

     1 1 1 1 . 1 Ω BB AB BB AA AB AA                          Since 0 Ω1
 and

BB AA

  , we can conclude than 0 1.

 

Second, we want to show that 1 



Let's prove it by contradiction. Then we make the hypothesis that 1   .

So we must have:

         11 1 11 1 1 1 1 (1 ) BB AB BB AB AA AB AA AB                                                    1 1 1 1 1 (1 )
. We can then rewrite the inequality: Appendix 4. Impact of the weights on the optimal value of industry share We want to assess the impact of the weights in the total utility function on the optimal value of the industry share. To do this, we normalize the weight for region A to 1 and give a weight  to region B. This way, the indirect utility function to be maximized is then:

             1 1 1 1 (1
  1 AB W V V                   max max 1 1 1 AA AB AB BB W W                                   
The first order condition gives us the following value: .

We want to compare Eq  with the optimal share of firms o  that can be rewritten as: To do so, we make the hypothesis that:

()

BB BB BB AB AA AB BB AB              
This would mean that :

    AA AB BB AB BB AB                      1 0 AA AB BB AB                  1 ( )

  The value of k then given by the following expression:

  have been calibrated, we must set the values of other parameters, that are not always known (for instance  ). Others, like the interregional transportation cost, can be found in UNCTAD(2004) .

Table 2 . Other parameters to run the simulation Other parameters

 2

	w wage	1
	 (elasticity of substitution)	6
	Share of transport cost in the price	22%
	of goods sold in the other region
	Iceberg cost AB 	1.28
		0.29
	AB
		0.6

Table 3 . Spatial equilibrium and spatial optimum Spatial equilibrium and optimum

 3

	Spatial equilibrium Eq 	0.75
	Intermediate parameter 	1.59
	Optimal concentration o 	0.69

Table 4 . Optimal taxes and subsidies on use of infrastructure to reach the optimum and their impact on the transportation costs Road pricing to reach the equilibrium

 4

	Optimal AA 	0.74
	Optimal tax	t	0.16
			AA
	Optimal			0.20
		AB
	Optimal tax	t	0.09
			AB
	Optimal			0.91
		BB
	Optimal subsidy	s	0.12
				BB

Table 5 . Impact of taxes on transportation costs Share of the transport cost in the price of a shipped product Before tax After tax

 5

	Interregional	22%	27%
	Intraregional in the center	2%	6%

  And this is false. So by contradiction, we can say that 1  . Since we have shown that  

										 θ
	Third, with 11    	, then we have		o 	[0,1]
	We know that 1 					1	  		, we now have:
										BB	AB
				.
	BB		AB
	This allows us to say that:
						λ	o	          BB AB BB AB AA AB ζ ζζ	[0,1].
				)  1 AA AB BB AB                   	AA BB        AB AB	1  	1	AA     	AB    BB 	AB	.



  First, we can compare the second parts of the two equations. Since   0

		o	( BB     )      				( AB     			)
			AA	AB	BB	AB	AA	AB	BB	AB
											 	  and
											BB	AB
	0
			 	AB	AB         AB    			.
			AA	AB		AA		AB		BB	AB
	Second, we want to compare the first parts of the equations. Let's prove by contradiction that
				BB   					() BB      
			BB	AB		AA		AB		BB	AB

AB

  , then we can observe that:

()

  This is a contradiction since we know that 01   . So this means that we have:These two inequalities lead us to the conclusion that: 

		BB    BB	AB			AA			AB	() BB BB AB      	.
										Eq	o
		AA      AB    	
								BB	AB
							1. 



λλ .

2 André de Palma and Stef Proost acknowledge the support of the French PREDIT research programme n° 09 MT CV 14 and of the Sustaincity project (FP7-EU)

Appendix 5. Proof of Proposition 3

The spatial equilibrium is given by the value Eq  :

      


P Hamon

M Gautier

P Garrec

A Janot

RAM

Keywords: robot, modeling, identification, friction I

Usually, the joint transmission friction model for robots is composed of a viscous friction force and of a constant dry friction force. However, according to the Coulomb law, the dry friction force depends linearly on the load driven by the transmission. It follows that this effect must be taken into account for robots working with large variation of the payload or inertial and gravity forces, and actuated with transmissions as speed reducer, screw-nut or worm gear. This paper proposes a new inverse dynamic identification model for n degrees of freedom (dof) serial robot, where the dry friction force is a linear function of both the dynamic and the external forces. A new identification procedure groups all the joint data collected while the robot is tracking planned trajectories with different payloads to get a global least squares estimation, in one step, of inertial and new friction parameters. An experimental validation is carried out with a 1 dof prismatic joint composed of a Star high precision ball screw drive positioning unit, which allows large and easy variations of the inertial and gravity forces.

INTRODUCTION

The usual off-line dynamic identification of robots uses the inverse dynamic identification model (IDIM) which is linear in relation to the dynamic parameters, and uses least squares (LS) technique. This procedure has been successfully applied to identify inertial and friction parameters of a lot of prototypes and industrial robots [START_REF] Atkeson | Estimation of Inertial Parameters of Manipulator Loads and Links[END_REF]- [START_REF] Kawasaki | Terminal-link parameter estimation and trajectory control of robotic manipulators[END_REF]. An approximation of the kinematic Coulomb friction, ( ) C F sign q ɺ , is widely used to model friction force at non zero velocity q ɺ , assuming that the friction force C F is a constant parameter. It is identified by moving the robot without load (or external force) or with constant given payloads [START_REF] Khalil | Identification of the payload inertial parameters of industrial manipulators[END_REF]. However, the Coulomb law suggests that C F depends on the transmission force driven in the mechanism. It depends on the dynamic and on the external forces applied through the joint drive chain. Consequently for robots with varying load, the identified IDIM is no more accurate when the transmission uses industrial speed reducer, screw-nut or worm gear because their efficiency significantly varies with the transmitted force. The significant dependence on load has been often observed for transmission elements [START_REF] Gogoussis | Coulomb friction effects on the dynamics of bearings and transmissions in precision robot mechanisms[END_REF]- [START_REF] Garrec | ABLE, an Innovative Transparent Exoskeleton for the Upper-Limb[END_REF] through direct measurement procedures. Moreover, the mechanism efficiency depends on the sense of power transfer leading to two distinct sets of friction parameters. This paper presents a new inverse dynamic identification model for n degrees of freedom (dof) serial robot, where the dry friction force C F is a linear function of both the dynamic and the external forces, and with asymmetrical behavior depending on the signs of the joint force and velocity. A new identification procedure is proposed. All the joint position and joint force data collected in several experiments, while the robot is tracking planned trajectories with different payloads, are concatenated to calculate a global least squares estimation in one step of both the inertial and the new friction parameters.

An experimental validation is carried out on a ball screw drive prismatic joint. This simple study case allows large and easy variations of the gravity and inertial forces. It is very suitable to study the effect of large variation of the load driven by the transmission on the friction force. The results obtained with the usual and the new friction model are compared.

II. USUAL INVERSE DYNAMIC MODELING AND IDENTIFICATION

A. Modeling

In the following, all mechanical variables are given in SI units in the joint space. All forces, positions, velocities and accelerations have a conventional positive sign in the same direction. That defines a motor convention for the mechanical behavior.

The dynamic model of a rigid robot composed of n moving links is written as follows [START_REF] Khalil | Modeling, identification and control of robots[END_REF]:

dyn in f ext = + + τ τ τ τ (1)
where: • dyn τ is the (nx1) vector of dynamic forces due to the inertial, centrifugal, Coriolis, and gravitational effects:

( ) ( ) ( ) dyn , = + + τ M q q C q q q Q q ɺɺ ɺ ɺ (2)
where q, q ɺ and q ɺɺ are respectively the (nx1) vectors of generalized joint positions, velocities and accelerations, M(q) is the (nxn) robot inertia matrix, ( , ) C q q ɺ is the (nxn) matrix of centrifugal and Coriolis effects, Q(q) is the (nx1) vector of gravitational forces.

• in τ is the (nx1) input torque vector of the drive chain: ( )

0 in f f f = - τ g v v (3)
where v f is the (nx1) vector of current references of the current amplifiers, 0 f v is a (nx1) vector of amplifiers offsets, g f is the (nxn) matrix of the drive gains,

f i t = g NG K (4)
N is the (nxn) gear ratios matrix of the joint drive chains ( m = q Nq ɺ ɺ with m q ɺ the (nx1) velocities vector on the motor side), G i is the (nxn) static gains diagonal matrix of the current amplifiers, K t is the (nxn) diagonal matrix of the electromagnetic motor torque constants.

• f τ is the (nx1) vector of the loss force due to frictions. Usually, it is approximated with a viscous friction and a dry friction:

( )

f V C Coff = - - - τ F q F sign q F ɺ ɺ (5)
where F V is the (nxn) diagonal matrix of viscous parameters, F C is the (nxn) diagonal matrix of dry friction parameters, and sign(.) denotes the sign function, F Coff is a (nx1) vector of asymmetrical Coulomb friction force between positive and negative velocities. This friction model is linear in relation to F V and F C (Fig. 1.a). • ext τ is the (nx1) external forces vector in the joint space.

Thus (1) becomes ( 6):

( ) ( )

( ) 0 dyn ext f f V C Coff f f out V C off - = - - - + ⇔ = - - - τ τ g v F q F sign q F g v τ τ F q F sign q τ ɺ ɺ ɺ ɺ (6)
where out dyn ext = τ τ τ is the output force (load force) of the drive chain,

0 off Coff f f = + τ F g v is an offset force that regroups
the amplifier offset and the asymmetrical Coulomb friction coefficient, and

f f = τ g v is the motor force without offset.
Then (1) can be rewritten as the inverse dynamic model (IDM) which calculates the motor torque vector τ as a function of the generalized coordinates:

( ) ( ) ( ) ( ) ( ) C V off ext out C V off , = + + + + + - = + + +
τ M q q C q q q Q q F sign q F q τ τ τ F sign q F q τ ɺɺ ɺ ɺ ɺ ɺ ɺ ɺ (7)

B. Identification

The choice of the modified Denavit and Hartenberg frames attached to each link allows to obtain a dynamic model linear in relation to a set of standard dynamic parameters St χ [START_REF] Gautier | Identification of robot dynamics[END_REF], [START_REF] Khalil | Modeling, identification and control of robots[END_REF]:

( ) St St , , = τ D q q q χ ɺ ɺɺ (8)
where ( )

St

, , D q q q ɺ ɺɺ is the regressor and St χ is the vector of the standard parameters which are the coefficients XX j , XY j , XZ j , YY j , YZ j , ZZ j of the inertia tensor of link j denoted j J j , the mass of the link j called m j , the first moments vector of link j around the origin of frame j denoted j M j = [MX j MY j MZ j ] T , the friction coefficients F Vj , F Cj , the actuator inertia called Ia j , and the offset τ offj . The velocities and accelerations are calculated using well tuned band pass filtering of the joint position [START_REF] Gautier | Identification of robot dynamics[END_REF].

The base parameters are the minimum number of parameters from which the dynamic model can be calculated. They are obtained by eliminating and by regrouping some standard inertial parameters [START_REF] Gautier | Direct calculation of minimum set of inertial parameters of serial robots[END_REF], [START_REF] Gautier | Numerical calculation of the base inertial parameters[END_REF]. The minimal inverse dynamic model can be written as:

( ) , , = τ D q q q χ ɺ ɺɺ (9)
where ( ) , , D q q q ɺ ɺɺ is the minimal regressor and χ is the vector of the base parameters.

The inverse dynamic model ( 7) is sampled while the robot is tracking a trajectory to get an over-determined linear system such that [START_REF] Gautier | Identification of robot dynamics[END_REF]:

( ) ( ) , , = + Y τ W q q q χ ρ ɺ ɺɺ (10)
with Y(τ) the measurements vector, W the observation matrix and ρ the vector of errors.

The LS solution χ minimizes the 2-norm of the vector of errors ρ. W is a (r×b) full rank and well conditioned matrix where e r N x n = , with N e the number of samples on the trajectories. The LS solution χ is given by: ( )

1 T T ˆ- + = = χ W W W Y W Y (11)
It is calculated using the QR factorization of W. Standard deviations i χ σ are estimated using classical and simple results from statistics. The matrix W is supposed to be deterministic, and ρ, a zero-mean additive independent noise, with a standard deviation such as:

( ) T E 2 r ρρ ρ σ = = C ρρ I ( 12
)
where E is the expectation operator and I r , the (r×r) identity matrix. An unbiased estimation of σ ρ is:

( ) 2 2 ˆr b ρ σ = - - Y Wχ ( 13
)
The covariance matrix of the standard deviation is calculated as follows:

T 2 T 1 χχ ρ E ( )( ) σ ( ) ˆˆˆˆ-  = - - =   C χ χ χ χ W W ( 14
) i 2 ˆˆˆii C χ χχ σ = is the i th diagonal coefficient of ˆχχ C . The relative standard deviation ri % χ
σ is given by:

ri i ˆˆi % 100 χ χ σ σ χ = (15)
However, experimental data are corrupted by noise and error modeling and W is not deterministic. This problem can be solved by filtering the measurement vector Y and the columns of the observation matrix W as described in [START_REF] Gautier | Dynamic identification of robots with power model[END_REF], [START_REF] Gautier | Extended Kalman filtering and weighted least squares dynamic identification of robot[END_REF].

III. NEW FRICTION MODELING AND IDENTIFICATION

In this section, we introduce a friction model dependent on the load, that is out τ .

A. Load-Dependent Friction Model

The Coulomb friction is still written ( ) C sign F q ɺ , with F C a (nxn) diagonal matrix. But here, for each link j,

( ) C j , j

F

(the (j,j) th element of the matrix C F ) depends linearly on the absolute value of the load of joint j which is out j τ (Fig. 1.b), [START_REF] Gogoussis | Coulomb friction effects on the dynamics of bearings and transmissions in precision robot mechanisms[END_REF]- [START_REF] Garrec | ABLE, an Innovative Transparent Exoskeleton for the Upper-Limb[END_REF]. As one can see in II.B, out j τ is a function of , , q q q ɺ ɺɺ and is linear in relation to base parameters.

Then the inverse dynamic model for each link j becomes:

( )

( ) ( ) j out j j out j j j V j , j j off j sign q q τ τ α τ β τ = + + + + F ɺ ɺ (16)
where the parameters j α and j β are constants to be identified.

These new parameters depend on the mechanical structure of the reducers used to actuate the robot.

For ease of understanding, the subscript j is omitted for all variables in the following to simplify the notation. The inverse dynamic model can be written as follows:

( ) ( )

out out V off sign q sign q F q τ τ α τ β τ = + + + + ɺ ɺ ɺ (17)
And with ( )

out out out sign τ τ τ = , one obtains: ( ) ( ) ( ) out out out V off sign sign q sign q F q τ τ ατ τ β τ = + + + + ɺ ɺ ɺ (18)
Thus, the IDM depends on the signs of out τ and q ɺ . With ( ) ( ) ( ) out out sign sign q sign q τ τ = ɺ ɺ , one defines 4 quadrants which can be grouped two by two (Fig. 2.a), depending on the sign of the output power denoted out out P q τ = ɺ . In the quadrants 1 and 3, P out is positive and the actuator has a motor behavior. In the quadrants 2 and 4, P out is negative and the actuator has a generator behavior which may save the power to the power supply, assuming a 4 quadrants amplifier. This model is valid for symmetrical friction. Generally, the friction is asymmetrical and α and β can take different values depending on the quadrant where the joint runs.

B. Power Sign-Dependent Friction Model

We present here 3 ways of modeling.

• 4 models (4x1Q) for 4 different quadrants

In the general case, the friction parameters α , β , and off τ depend on the signs of τ load and q ɺ in the frame ( q ɺ , load τ ), Fig. ( )

( ) ( ) ( ) load 1 out 1 V load 2 out 2 V load 3 out 3 V load 4 out 4 V 0 & q 0 1 F q 0 & q 0 1 F q 0 & q 0 1 F q 0 & q 0 1 F q τ τ α τ β τ τ α τ β τ τ α τ β τ τ α τ β > > ⇒ = + + +   > < ⇒ = - -+   < < ⇒ = + -+   < > ⇒ = - + +  ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ (19)
• 2 models (2x2Q) for quadrants identical two by two

In some cases, the friction is symmetrical with respect to the velocity and the 4 models can be simplified to 2 models.

out m out m V off m out g out g V off g P 0 1 sign q F q P 0 1 sign q F q τ α τ β τ τ α τ β τ > ⇒ = + + + +    < ⇒ = - + + +   ɺ ɺ ɺ ɺ (20)
0 q 0 < ɺ 0 q 0 < ɺ 0 q 0 > ɺ 0 q0 > ɺ b) q ɺ β -β out τ increases Friction a) q ɺ F c -F c Friction • 1 model (1x4Q) for 4 identical quadrants
In the case of symmetrical friction with respect to out τ and q ɺ , the models simplify to 1 model for 4 quadrants.

( ) ( )

out out out V off sign P sign q F q τ τ ατ β τ = + + + + ɺ ɺ (21)
Each modeling takes the load-dependency of friction into account. Starting from the (4x1Q) models, a simplification procedure based on the identification results leads to the simplest model. However, when the joint torque is low, the friction model should be more complex because of Coulomb friction resulting from internal preload and hysteresis. Then, one cannot use the quadrants. To simplify and connect the quadrants models, a relevant approximation consists in extending them in this area which can be considered as uncertain for the experimental identification. This simplification is illustrated in a ( τ , out τ ) graph (Fig. 2.b). It should be noticed that within this area, the mechanism is no longer transmitting power which is totally dissipated in friction losses.

C. Friction Identification Method

The friction models depend on the sign of out τ which is unknown. To overcome this problem, the samples of τ measurements are selected outside of the uncertain area (Fig. The (1x4Q) modeling can be written:

( ) ( ) out out V off sign P sign q F q τ τ ατ β τ = + + + + ɺ ɺ (23)
For the (2x2Q) modeling, 2 variables are introduced, P + and P -, defined by: ( )

1

IV. PROTOTYPE TO BE IDENTIFIED

The EMPS is a high-precision linear Electro-Mechanical Positioning System. Its main components are a Maxon DC motor which is current controlled by a four quadrant PWM amplifier, a Star high-precision low-friction ball screw drive positioning unit, and an incremental encoder. The backlash free ball screw drive is the gear converting the rotary motion of the motor to the linear carriage joint displacement. The EMPS is a standard configuration of a drive system for prismatic joint of robots, machine tools, haptic device… It is connected to a dSPACE digital control system for easy control and data acquisition using Matlab and Simulink software [22], [23]. In order to make easy variation of the gravity load, the EMPS can be fixed in vertical position alternatively with the z joint positive linear displacement in the gravity direction or in its opposite.

The inverse dynamic model is given by:

( ) ( ) ( ) ( ) ( ) C V off a C V off M q q Q q F sign q F q I mq mg F sign q F q τ τ τ = + + + + = + - + + + ɺɺ ɺ ɺ ɺɺ ɺ ɺ ( 27
)
where:

• I a is the inertia moment of all rotary elements in the drive chain (rotor of motor and encoder, ball, coupling units), • m is the mass of all translation moving elements (screw, carriage, payload), • g is the projection of the gravity acceleration on the z prismatic joint axis. • g = 9.81, when z axis stays in vertical position, oriented positive to the earth, • g = -9.81, when z axis stays in vertical position, oriented positive to the sky.

All variables and parameters are given in SI units on the joint linear displacement side (carriage side).

The identification process has been performed for the four different cases: first with the usual model where F C is constant, then with the 3 friction models depending on the load. We have observed that the modeling with four identical quadrants (1x4Q) was insufficient because the friction is asymmetrical here. Moreover, the modeling with four different quadrants (4x1Q) was not indispensable as the parameters are very close for the two motor quadrants on one hand and for the two generator quadrants on the other hand. That is the reason why only the modeling with quadrants identical two by two (2x2Q) is detailed here and compared with the usual model. The method for the others is very similar.

V. EXPERIMENTAL IDENTIFICATION

A. Data Acquisition

The identification of dynamic parameters is carried out with and without payloads: five different additional masses can be fixed on the carriage. To excite properly the friction parameters to be identified, triangular trajectories were used. The sample acquisition frequency for joint position and current reference (drive force) is 5 KHz.

The estimation of q ɺ and q ɺɺ are carried out with pass band filtering of q consisting of a low pass Butterworth filter and a central derivative algorithm. The Matlab function filtfilt is used. This is a zero-phase forward and reverse digital filtering. We calculate the drive force using the relation:

f f g v τ = (28)
where v f is the current reference of the amplifier current loop, and g f is the gain of the joint drive chain, which is taken as a constant in the frequency range of the robot because of the large bandwidth (700 Hz) of the current loop [START_REF] Restrepo | Calibration of drive chain of robot joints[END_REF].

In order to cancel high frequency ripple in v f (and τ ) , the vector Y and the columns of the observation matrix W are both low pass filtered and decimated. This parallel filtering procedure can be carried out with the Matlab decimate function [START_REF] Swevers | Optimal Robot excitation and Identification[END_REF], [START_REF] Kawasaki | Terminal-link parameter estimation and trajectory control of robotic manipulators[END_REF].

B. Identification

To identify the load-dependant friction, measurements with known payloads must be used and one needs the relation:

0 a m m m = + ( 29
)
where:

• m 0 is the unknown mass of the translational elements, with the carriage free of additional mass, • m a is one of the 6 additional masses, fixed on the carriage, with accurate weighted values: 0 kg, 1.05 kg, 3.0266 kg, 4.7882 kg, 9.9162 kg, and 14.704 kg.

Hence for the samples ( )

k τ with the additional mass m a(k) , (24) becomes: ( ) ( ) ( ) ( ) ( )
k a 0 0 a k C V off I m q m g m q g F sign q F q τ τ = + - + -+ + + ɺɺ ɺɺ ɺ ɺ (30)
At a first step, we proceed separately for each (k) experiment to get 6 different usual identifications. Keeping only the samples at average constant velocities without acceleration, the load corresponds to the gravity effect. One uses usual LS method as described in II.B. The variation of F C as a function of the 6 payload values is given in Fig. 6 (see the asterisk*), in order to show the linear dependency on load.

At a second step, to identify the usual model with all samples, one distinguishes the weighed mass m aw and the mass m ae estimated by the identification. Thus, the usual model is written:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ae k k a 0 0 aw k aw k C V off m I m q m g m q g ... m ... F sign q F q τ τ = + - + -+ + + ɺɺ ɺɺ ɺ ɺ (31)
Then, the sampled measurements, for k from 1 to 6, are concatenated using the m aw(k) weighed value corresponding to each experiment (k), to get the linear system:

usual usual usual usual = + Y W χ ρ ( 32
)
with the measurements vector, the observation matrix, and the vector of base parameters below:

T T T usual (1) (6)
...

  = =   Y τ τ τ ( 33
) usual aw ( ) ( ) = - -     W q g m q g sign q q 1 ɺɺ ɺɺ ɺ ɺ (34)
T

ae usual a 0 0 C V off aw m I m m F F m τ   = +     χ ( 35
)
At a third step, the proposed model is identified with:

new new new new = + Y W χ ρ (36)
with the measurements vector, the observation matrix, and the vector of base parameters defined as follows:

new usual = = Y Y τ (37) new aw aw ( ) ( ) ( ) ( )
+ + + + + - - - - -  = - -   - - 
W P q P g P m q g P sign q P ... P q P g P m q g P sign q P q ɺɺ ɺɺ ɺ ɺɺ ɺɺ ɺ ɺ

(38) T ( )( ) ( ) ( ) ( )( ) ( ) ( ) new m a 0 m 0 m m off m g a 0 g 0 g g off g V 1 I m 1 m 1 ... 1 I m 1 m 1 F α α α β τ α α α β τ  = + + + +   - + - -  χ ( 39
)
Here + P and - P are diagonal matrices, with: ( ) ( )

i i i ,i i ,i 1 sign 1 sign , 2 2 + - + - = = P P P P ( 40
)
The two models are compared using exactly the same identification method with the same measurements.

C. Results

The values identified with usual IDM and OLS regressions are given in Table I and those with the new IDM in Table II. For each model, Fig. 4 and Fig. 5 present a direct validation comparing the actual τ with its predicted value Ŵ χ . In Fig. 6, it can be seen that the variation of the usual F C is the mean between motor and generator values, except for the 2 first masses because the load is to low and the joint works only in motor mode (see the uncertain area Fig. Table III shows that the new model does not improve the residual if the load variation is too small. Indeed, in this case (additional mass less than 3 kg), the joint runs mostly in motor behaviour and the usual model is then sufficient. For additional mass greater than 5 kg, the load variation is large enough to justify the use of the proposed model: one observes a decrease of 30% in the relative norm of errors.

VI. DISCUSSION

The new IDM with a load-dependent joint friction model brings a substantial improvement for joint whose load can vary significantly. Robots carrying important masses or with large variation of inertial and gravity forces are considered. Moreover, the identification process is the same for both usual and new models, because the new model remains linear in relation to the parameters.

The main difficulty is to distinguish the different behaviours, motor and generator, but a solution has been proposed along this paper. However, the measurements have to be more exciting than usual. Each test has to be done with different loads to highlight the effect on the friction variations. So, this identification protocol is more time-consuming and the setting up must be adapted for the measurements with additional masses. Moreover, for robots with small load variation, the joint actuates only in one quadrant: the parameters m α and g α of the new model are not excited enough and the results are not better than with the usual. Then, this one will be accurate enough for low load variation.

One has to consider the load variation rate and the type of the transmission to choose the appropriate model, before starting modeling and identification.

VII. CONCLUSION

In this paper, a new model of friction with a loaddependence has been presented. The identification process has been adapted to the new inverse dynamic identification model. This method was validated experimentally on a 1 dof prismatic joint robot. In addition, we have established that this technique can be applied to multi dof robots because we get a linear system with respect to friction and dynamic parameters. The proposed model is accurate and required for robots with large load variation.

Future works concern the application of this IDIM to identify the parameters of a multi dof robot. Several types of transmission will be studied, as well as a possible velocitydependence of dry friction. Once an accurate model will be identified for multi dof robots, it will be used for torques monitoring and collision detection.

Figure 1 .

 1 Figure 1. a) Usual friction model with constant FC. b) Parametric effect of the load on friction model.

Figure 2 .

 2 Figure 2. a) Four quadrants frame ( ) out q , τ ɺ

  e r a t o r m o to r g e n e r a t o r

  2.b) in order to get the same sign for out τ and τ . This allows to get the sign of P out with:

Figure 3 .

 3 Figure 3. EMPS prototype to be identified.

Figure 4 .

 4 Figure 4. Direct validation performed with usual IDM

Figure 6 .

 6 Figure 6. Evolution of Coulomb friction FC and comparison between the two models.

TABLE I

 I

					.	IDENTIFIED VALUES WITH USUAL IDM
	Parameters		Identified Values		Standard deviation * 2		Relative deviation
		Ia + m0			64.799			0.477			0.4 %
		m0			1.047			0.012			0.6 %
		mae/maw		1.025			0.002			0.1 %
		FC			38.277			0.237			0.3 %
		FV			396.550			2.894			0.4 %
		τoff			-7.935			0.078			0.5 %
			TABLE II.		IDENTIFIED VALUES WITH NEW IDM
	Parameters		Identified Values		Standard deviation * 2		Relative deviation
	(1+αm)(Ia+m0)		65.850			0.375			0.3 %
		(1+αm)m0		0.821			0.008			0.5 %
		(1+αm)			1.174			0.001			0.1 %
		βm			31.337			0.171			0.3 %
		τoffm			-8.415			0.059			0.4 %
	(1-αg)(Ia+m0)		68.339			0.692			0.5 %
		(1-αg)m0		0.711			0.182			12.8 %
		(1-αg)			0.831			0.004			0.2 %
		βg			21.780			1.803			4.1 %
		τoffg			-7.070			0.122			0.9 %
		FV			409.285			1.983			0.2 %
		400
		300
		200
		100
	Input τ (N)	0
		-100
		-200
		-300										Measurement Estimation
												Error
		-400	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
								Samples				x 10 4

  In direct validation, it is shown that the predicted torque is improved with the new IDM. It can be seen that, for the EMPS, TableIIIpresents the relative norm of errors ρ Y for the two models and for two sets of experiments.

	also why the parameter ( 1	-	α	g	)	0 m	is not excited enough and
	not well identified.
		400
		300
		200
		100
	Input τ (N)	0
		-100
		-200
		-300											Measurement Estimation
													Error
		-400	0	0.2	0.4	0.6	0.8	1		1.2	1.4	1.6	1.8	2
									Samples	x 10 4
				Figure 5. Direct validation performed with new IDM.
	m α and g α are very close whereas m β and g β are
	significantly different.
	Moreover,
	2.b). That explains

TABLE III

 III

	.	RELATIVE NORM OF ERRORS WITH BOTH MODELS
	Measurements used for the identification	Relative norm of errors with the usual model	Relative norm of errors with the new model
	Samples without additional mass
	and samples with the additional	0,0916	0,0908
	masses of 1,05 kg and 3,0266 kg
	All samples (without and with each additional mass)	0,09996	0,0679