L Nehaoua
email: nehaoua@ibisc.univ-evry.fr

L Nouvellière
email: nouvelliere@ibisc.univ-evry.fr

S Mammar
email: mammar@ibisc.univ-evry.fr

Dynamics modeling of a Two-wheeled vehicle using Jourdain's principle

Keywords:

This paper describes a modeling technique for developing the motion equation of a motorcycles vehicle. Based on the Jourdain's principle approach, the derived model presents a nice analytical formulation with recursive Jacobian matrices computation. The synthesized model takes in consideration the main forces and moments affecting the behavior of the motorcycle and allowing the simulation of 11DOF.

I. INTRODUCTION

Riding a two-wheeled vehicle is not simple and intuitive like driving a car vehicle. Involved accelerations, critical equilibrium, aerodynamic effects, trajectory anticipation are, among others, many points that make driving a motorbike to be a complicated exercise.

To predict the behavior of such a system, one method is to model it. In the study of mechanical systems, modeling usually leads to a set of differential algebraic equations (DAE) derived from laws of mechanics. In this field, two approaches are proposed: the analytical approach and the numerical one. In this later, many tools were developed such MSC ADAMS [START_REF]ADvAnced Multi-body Simulation[END_REF], simMechanics [2], Dymola [START_REF] Modeling | [END_REF], etc.. The alternative approach is the development of symbolic equations of motion. It is very time consuming method, however, it offers more flexibilities in the system description. In this case, using a dedicated packages implemented in a computer algebra software, like TsiProPac [4] or MBsymba [START_REF] Lot | A Symbolic Approach for Automatic Generation of the Equations of Motion of Multibody Systems[END_REF], is highly recommended.

Many works have been addressed the modeling of the twowheeled vehicles. Indeed, the first study on the stability and balance of a bicycle vehicle were performed by Whipple [START_REF] Whipple | Stability of the motion of a bicycle[END_REF]. This bicycle is represented by two bodies linked via the steering mechanism. Wheels are modeled as circular thin disks where the gyroscopic effect is taken into account. After Whipple, other work has been undertaken [START_REF] Bower | Steering and stability of single track vehicles[END_REF], [START_REF] Pearsall | The stability of the bicycle[END_REF], [START_REF] Kondo | Theoretical study on the ruuning stability of twowheelers[END_REF]. In all these studies, the tire-road interaction has been neglected. In 1971, Sharp integrated tires' efforts [START_REF] Sharp | The stability and control of motorcycles[END_REF]. The objective of this study is to examine the effect on stability of a variation in the vehicle's parameters. His model is represented as a set of two rigid bodies connected by the steering mechanism, which can simulate 4 Degrees of Freedom (DOF).

So far, in all these models, the longitudinal velocity was assumed to be constant and the vehicle suspension movements were ignored. It was interesting to examine how these factors influenced the frequency characteristics of the different stability modes. These aspects have been investigated in a further work of Sharp in 1974 [START_REF] Sharp | The influence of frame flexibility on the lateral stability of motorcycles[END_REF]. In 1978, Weir made the first approach to stabilize a motorcycle by including a human in the simulation loop [START_REF] Weir | Motorcycle handling[END_REF]. In 1994 Sharp published a paper, in which, the motorcycle is represented by a set of four rigid bodies to simulate 8DOF. The originality of this model lies in the integration of flexibilities and the rider tilting motion [START_REF] Sharp | Vibrational modes of motorcycles and their design parameter sensitivities[END_REF].

In the same vein, a model incorporating the suspension travels and a more complex representation of the tire-road interaction is presented in [START_REF] Sharp | A Motorcycle Model for Stability and Control Analysis[END_REF]. The advances made by this model have allowed to explore the dynamics of such a vehicle for a large motion around equilibrium. In this case, the interaction between the in-plane and the out-of-plane motion is best described. More recently, the project FastBike presents a nonlinear model of a motorcycle described as a set of 6 bodies and simulating 11DOF. Each body is identified by its natural coordinates and the assembly is done by introducing the geometric constraints [START_REF] Cossalter | A motorcycle multibody model for real time simulation based on the natural coordinates approach[END_REF]. Finally, a modeling technique based on the recursive Newton-Euler approach is adapted to derive a motorcycle dynamics model intended for a two-wheeled simulator application [START_REF] Hima | Suitable Two Wheeled Vehicle Dynamics Synthesis for Interactive Motorcycle Simulator[END_REF].

In the present paper, a modeling technique based on the Jourdain's principle is used to develop the motion equation of a motorcycles vehicle.

II. JOURDAIN'S PRINCIPLE

A. Preliminaries

Consider a mechanical system of n bodies described by m natural coordinates. We assume that the system dynamics will be expressed into a moving reference frame ℜ v . The configuration of the moving reference ℜ v with respect to (w.r.t) the inertial reference ℜ o is given by the position vector rov = [w, y, z] T and orientation vector qv = [ϕ, θ, ψ] T . The position of a body i w.r.t the inertial reference ℜ o expressed in this latter is given by:

ro oi = ro ov + ro vi (1)
and:

ωo ov = R o→v ωv ov → ǭv ov = ωv ov (7)
By projecting equations ( 3) and ( 5) into the ℜ v reference and by using ( 6) and [START_REF] Bower | Steering and stability of single track vehicles[END_REF] we can define the linear and angular acceleration of the each body i as following:

āv oi = vv ov + ωv ov × rv vi + rv vi + ωv ov × ( ṙv vi + vv oi ) (8) and:
ǭv oi = ωv ov + ωv ov × ωv vi + ωv vi (9)
Equation ( 8) and ( 9) are useful to deduce Jacobian matrices, as it will be demonstrated in the next sub-section.

B. Jourdain's principle

Suppose that the kinematics of the mechanical system is characterized by k generalized speeds θ = [v v ov , ωv ov , qs ] corresponding to m natural coordinates χ = [r o ov , qv , qs ]. Consequently, the linear and angular velocity of each body can be expressed as followings (henceforth, all vectors are expressed in ℜ v reference unless otherwise stated):

voi = voi ( θ, χ) ωoi = ωoi ( θ, χ)
Then, the linear acceleration is expressed by:

āoi = ∂v oi ( θ, χ) ∂ θ θ + ∂v oi ( θ, χ) ∂χ χ
and the same equation for ωoi can be deduced. The partial derivatives are called Jacobian matrices. Let the position and the angular rates of the i-th body w.r.t the ℜ v to be:

ṙv vi = ᾱi qs ωv vi = βi qs (10)
Now, it is easy to construct the i-th Jacobian matrix as shown in the following two equations:

∂v oi ∂ θ = I 3 -r vi ᾱi (11)
∂ ωoi ∂ θ = 0 3 I 3 βi
Now, all kinematics quantities are defined, then we introduce the principle of the virtual power called also Jourdain's principle. Jourdain's principle states that the virtual power done by motion compatible constraint efforts Fi,c / Mi,c is null [START_REF] Rill | Simulation von Kraft-fahrzeugen[END_REF], then:

n δv T oi Fi,c + δ ωT oi Mi,c = 0 (12)
Next, the Newton/Euler dynamics principle allows to express the constraint efforts w.r.t to the external applied efforts Fi,e / Mi,e as:

m i āoi = Fi,c + Fi,e I i ǭoi + ωoi × I i ωoi = Mi,c + Mi,e
By replacing this equations into [START_REF] Weir | Motorcycle handling[END_REF], and after some algebraic manipulations, we deduce the dynamics motion expression as following:

M θ = Q ( 13
)
where M is the mass matrix obtained by the direct calculation of Jacobian matrices:

M = m i ∂v T oi ∂ θ ∂v oi ∂ θ + ∂ ωT oi ∂ θ I i ∂ ωoi ∂ θ (14)
and the vector Q of the generalized efforts can be split into two parts Q = Qe -Qr , Qe is the generalized efforts vector associated to the external applied efforts and Qr is the generalized efforts vector associated to the acceleration residual terms in acceleration equations [START_REF] Pearsall | The stability of the bicycle[END_REF][START_REF] Kondo | Theoretical study on the ruuning stability of twowheelers[END_REF]. Hence:

Qr = m i ∂v T oi ∂ θ ∂v oi ∂ χ χ (15)
+ ∂ ωT oi ∂ θ ∂ ωoi ∂ χ χ + ωoi × I i ωoi (16)
Qe = ∂v T oi ∂ θ Fi,e + ∂ ωT oi ∂ θ Mi,e (17)

III. VEHICLE KINEMATICS

In this section, a direct application of the Jourdain's principle (section II) to the motorcycle modeling is exposed. For this, the kinematics of each body w.r.t the moving reference ℜ v is demonstrated.

At first, we highlight the different motion reference frames used in the current modeling work. As shown in figure [START_REF]ADvAnced Multi-body Simulation[END_REF], the moving reference ℜ v is obtained after two successive rotations w.r.t the inertial reference ℜ O . First, a yaw rotation ψ around ℜ o z-axes is done to give an intermediate reference

ℜ ψ . Next, ℜ v is obtained by a roll rotation ϕ around ℜ ψ x- axes. The origin v of ℜ v is defined by χv = [r o ov , qv ].
The motorcycle is considered as a set of six bodies (Fig. 1). The main body G r , the front assembly G f including the handlebar, steering column and the higher fork part, the unsprung front assembly G u , the rear swing arm G s and the two wheels G Rf , G Rr . The mass of the front wheel is included in the unsprung front assembly, only its rotation will be considered as an independent DOF. The rider motion is not taken into account and is considered to be rigidly attached to the main body G r .

In order to completely define the configuration of the motorcycle, a set of 11DOF has been considered. A full characterization of motorcycle's kinematics includes :

• the longitudinal, lateral and vertical position (x, y, z) of the main body, • the roll, pitch and yaw orientation (ϕ, θ, ψ) of the main body, • the steering angle δ,

• the front suspension travel λ f and the swing arm rotation µ, • and finally, the front and rear tires spinning, ξ f and ξ r . In the following subsections, the kinematics of each body is presented.

A. Main body G r

The main body G r has 1DOF w.r.t ℜ v which represents the vertical displacement z. Then, its position vector is: [START_REF] Sharp | The stability and control of motorcycles[END_REF], we conclude that: ᾱGr = ēz 0 3,6 βGr = 0 3,7

rv vGr = 0 0 z + h r T → ṙv vGr = żē z From equation

B. Front body G f

The front body G f is belonging to the fixed-body reference frame ℜ b and has 2DOF w.r.t the main body G r . ℜ b frame is obtained by a three successive rotations namely, pitch θ about the y-axis of ℜ v , caster ǫ about the same axes and the steer δ about the ℜ b z-axis. The position vector of the front body w.r.t ℜ v is:

rv vG f = rv vGr + R v→g (r v GrB + R g→b rb BG f )
Knowing that vectors rb BG f and rv GrB are constant, we can write :

ṙv vG f = ṙv vGr + ωv vg × rv GrG f + ωv gb × rv BG f
where ωv vg = θē y and ωv gb = δē δ . In this case, ᾱG f = ᾱGr + ᾱ′ G f where:

ᾱ′ G f = 0 3,1 ēy × rv GrG f ēδ × rv BG f 0 3,4
In the other hand, the angular rates is ωv vG f = ωv vg + ωv gb , and hence:

βG f = 0 3,1 ēy ēδ 0 3,4
It can be seen that the advantage of this method, is that the computation of Jacobian matrices can be obtained recursively. At each step, only the contribution of the considered body is calculated.

C. Front unsprung body G u

The front unsprung body has 1DOF w.r.t G f representing the suspension travel λ f . Its position vector is:

rv vGu = rv vG f + R v→g R g→b rb G f Gu
Herein, rb G f Gu = (λ f + λ f 0 )ē z where λ f 0 is the initial fork suspension length. By differentiation of the position equation, we can write :

ṙv vGu = ṙv vG f + (ω v vg + ωv gb )r v G f Gu + λf ēδ
Consequently, ᾱGu = ᾱG f + ᾱ′ Gu , where:

ᾱ′ Gu = 0 3,1 ēy × rv G f Gu ēδ × rv G f Gu ēδ 0 3,4
Finally, since there is no rotation between the two front bodies: βGu = βG f .

D. Swing arm body G s

The swing arm body is belonging to the fixed-body reference frame ℜ p and has 2DOF w.r.t the main body G r . ℜ p frame is obtained by two successive rotations, pitch θ and swing arm pitch µ around the y-axis of ℜ v . Its position vector is expressed by::

rv vGs = rv vGr + R v→g (r g GrP + R g→p rp P Gs
) and the corresponding position rate is :

ṙv vGs = ṙv vGr + ωv vg × rv GrGs + ωv gp × rv P Gs
where ωv gp = μē y . Then, ᾱGs = ᾱGr + ᾱ′ Gs , where:

ᾱ′ Gs = 0 3,1 ēy × rv GrGs 0 3,2 ēy × rv P Gs 0 3,2
In the other hand, the angular rates is ωv vGs = ωv vg + ωv gp , and hence: βGs = 0 3,1 ēy 0 3,2 ēy 0 3,2

E. Front and rear wheels body

The two wheels has one independent generalized velocity w.r.t to their respective previous body (it means G u for the front wheel and G s for the rear one), which consist simply on tire rotation around wheel spin axes ēyR . It follows that ωv vGRr = ωv vGs + ξr ēy this leads to :

βGRr = βGs + 0 3,6 ēyR,r
The same reasoning can be followed for the front wheel:

βG Rf = βGu + 0 3,5 ēyR,f 0 3,5

At now, the whole motorcycle kinematics is defined. Matrix mass M and the first part of generalized efforts Qr (equations 14, 15) are fully defined. In the next section, we look for the development of the generalized efforts vector associated to external applied efforts.

IV. VECTOR OF GENERALIZED EFFORTS

To express generalized efforts vector associated to external applied efforts, one can make use of equation ( 16). First, it is essential two establish the set of all efforts exerted on each body. This method may be a tedious one especially one the number of applied efforts is much high. So, it is most convenient to find the virtual power done by each effort and find the associated generalized effort.

As an example, consider a body i under a force vector Fe applied at a point j ∈ i and a moment vector Fe around an axes belonging to this body. The resulting virtual power can be expressed as:

δP = δv T oj Fe + δ ωT i Me = δ θT Q
By replacing δv oj and δ ωi by their Jacobians, we obtain the following equation:

Q = ∂v oj ∂ θ T Fe + ∂ ωoi ∂ θ T Me (18)
The only difference w.r.t equation ( 16) is the use of the Jacobian related to the force point application velocity (v oj ) instead of that of the body center of mass voi . In the following, we use equation [START_REF] Rill | Simulation von Kraft-fahrzeugen[END_REF] for each external force/moment exerted on the motorcycle vehicle.

A. Front and rear suspension

The front suspension force Fsf is applied between the the front body G f and the front unsprung body G u . The contribution of this force in the total generalized effort is given by using equation [START_REF] Rill | Simulation von Kraft-fahrzeugen[END_REF], like following:

Qsf = ∂ ∂ θ voG f -voGu T Fsf (19)
By replacing each Jacobian by its expression and after some algebraic manipulations we get: Knowing that the force suspension vector can be written as Fsf = F sf ūsf , where, ūsf = ēδ is the unit vector, then:

Qsf =   0 3,3 (λ f + λ f 0 )ẽ T δ -ᾱ ′T   Fsf (20)
Qsf = -F sf ∂ λf ∂ϑ (21)
Doing the same for the rear suspension, we conclude :

Qsr =   0 8,1 (ē y × rv P s l ) T ūsr 0 2,1   F sr (22)
Equation ( 21) shows that the rear suspension force can be viewed as a stiffness torque about the swing arm pivot P .

B. Gravity forces

Each body is subject to the gravity force created by its own mass. The contribution of this force in Q is directly given by:

Q g = n m i ∂v oi ∂ϑ T ḡ (23)

C. Rider action

The rider's torque τ r exerted on the motorcycle's handlebar is considered to be the most important input for the vehicle driving. This torque acts around the steering axis and hence jointly on the main body G r and front body G f . The contribution of this torque in Q is defined by:

Qτr = ∂ ∂ θ ωG f -ωGr T τ r ēδ =   0 8,1 τ r 0 4,1   (24)

D. Propulsion and braking torques

For the front wheel, braking force acts between the wheel rim and unsprung assembly G u which contains braking disk. For the rear wheel, the braking torque is exerted between the rear wheel body and the swing arm assembly. In addition to the propulsive torque, the contribution of these torques is expressed by:

QDB = ∂ ∂ θ ωG Rf -ωGu T τ bf ēyR,f + ∂ ∂ θ (ω GRr -ωGs ) T (τ br + τ D )ē yR,r
After simplification we get:

Gu QDB =   0 11,1 τ bf τ br + τ D   (25)
In the next section, we expose the tire-road interaction to fully defining the vector of the generalized efforts Q.

V. TIRES KINEMATICS AND DYNAMICS

A. Tires kinematics

To describe the tire motion, a new reference frame ℜ T (C, ēxT , ēyT , ēzT ) is introduced at the contact point C. ēzT is the normal vector to the road surface. Vector ēxT is obtained as ēxT = ēyR × ēzT and ēyT completes the reference axis. As shown in figure 3, the camber angle is an important variable in the motorcycle safety studies. This variable is expressed by the following equation:

sin γ = ēT yR ēzT (26)
The position of the contact point w.r.t the inertial frame is expressed by:

ro oC = ro oGR -R rm ēo zR -R cr ēo zT (27)
Where, R rm is the wheel rim radius, R cr is the dynamic tire crown radius and ēzR = ēxT × ēyR . By differentiating this equation, one obtains the contact point velocity vector, as:

vo oC = vo oGR -R rm ėo zR -Ṙcr ēo zT (28)
Knowing that the contact point C is belonging always to the road surface (Fig. 4) then, R cr et Ṙcr are obtained by performing a product scalar of equations (26,27) by ēo zT . Once the tire's kinematics is established, the equivalent tire forces/moments wrench at the center of each wheel is given by the following expressions:

QT = j ∂v oGRj ∂ θ T FT + ∂ ωoGRj ∂ θ T MT (29)
With QT , the vector of the generalized efforts Qe is fully defined, and:

Qe = Qsf + Qsr + Qg + Qτr + QDB + QT (30)

VI. SIMULATION

To test the behavior of the developed motorcycle model, a linear PID controller is synthesized for the longitudinal speed regulation. This controller acts like an acceleration/braking torque directly applied on the rear wheel. The control torque is saturated by the maximum torque admissible by the engine at the rear tire axis. For this test, we drive the motorcycle from an initial speed of 10 km/h to a desired velocity of 72 km/h. Figure [START_REF] Lot | A Symbolic Approach for Automatic Generation of the Equations of Motion of Multibody Systems[END_REF] shows the control driving torque as computed by the PID controller and the longitudinal speed profile. Figure [START_REF] Whipple | Stability of the motion of a bicycle[END_REF] describes the resulting longitudinal acceleration and the front and rear suspension travels. Finally, figure [START_REF] Bower | Steering and stability of single track vehicles[END_REF] represents the vertical force F z applied at the contact point of each tire. the two hiccups in the speed vs time curve (and in the other figures) are due to the gear changes introduced by a simple logic switching (transmission model is not described in this paper). Although we have no experimental data, these simulations seems to be very correct in comparison with a previously published work in literature. The lateral dynamics simulation is very hard to carried out. Indeed, the motorcycle vehicle is highly unstable and a suitable controller must be synthesized.

VII. CONCLUSION

Based on the virtual power principle, called also Jourdain's principle, a motorcycle dynamic model is proposed. The development method is chosen for its simplicity and its fine analytical derivation with an acceptable calculation burden. A first test is presented to show the ability of the proposed model to simulate the longitudinal/vertical motion.Future works will focus on the lateral mode simulation by implementing a suitable controller especially for trajectory tracking like cornering, slalom and lane changes.

Fig. 1 .

 1 Fig. 1. The motorcycle geometrical description

Fig. 2 .

 2 Fig. 2. Suspensions configuration, (a) front, (b) rear

Fig. 3 .

 3 Fig. 3. Tire reference frame and camber

FT =F

  x ēxT + F y ēyT + F z ēzT MT = M x ēxT + M y ēyT + M z ēzT + FT × rCGR where F x = F x (κ) is the longitudinal force, F y = F y (ζ, γ) is the lateral force and F z = F z (δ z ) is the vertical force. M x is the tire torque about ēxT , M y is the rolling resistance torque and M z = M z (ζ, γ) is the alignment torque. This forces and torques are obtained from the longitudinal slip κ, the lateral slip ζ and δ z is the vertical tire deformation.

Fig. 4 .

 4 Fig. 4. Wheel and tire geometry: rim and crown

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Driving torque and the longitudinal speed

  VIII. ACKNOWLEDGMENT IX. APPENDIX• ē variables with over bar are vectors.• M calligraphic variables are matrices.• ēx = [1, 0, 0], ēy = [0, 1, 0], ēz = [0, 0, 1],• for an arbitrary rotation φ, we have:v→g = R θ,y , • R g→b = R ǫ,y R δ,z , • ēδ = R v→g R g→b e z , • ēv yR,f = R v→g R g→b ēy , • if x = [x, y, z] T , x =

This research was supported by funds from the NATIONAL RESEARCH AGENCY in PARTAGE project.


N Ait-Oufroukh

K Messaoudene

S Mammar

Dynamic Model of Steer-by-wire System for Driver Handwheel Feedback

Keywords: Steer-by-wire, Dynamic model, handwheel, tires friction

Recent advances towards steer-by-wire technology have promised significant improvements in vehicle handling performance and safety. This technology inspired from aeronautics, consists in replacing the traditional mechanical linkages and/or hydraulic power assist with electrical equivalents. This paper describes the adopted approach for development of a new dynamic model of a Steer-by-wire system with hydraulic assistance in the case of stopped vehicle. The aim is to establish dynamic equations of the forces and torques exchanged between handwheel and the front wheels during steering operations. Some simulation results are provided. The proposed model will be experimented using a prototype equipped with handwheel feedback system.

INTRODUCTION

The aim of electronics improvement is to enhance safety, reliability and driving comfort, while reducing the manufacturing costs. Advances m electronics have revolutionized many automotive engineering aspects, especially in the engine combustion management field and vehicle safety systems. The important benefit of applying electronic technology is to improve the performance of vehicle. However, only recently has the electronic begun to find its way into automotive steering systems in the form of electronically controlled variable assistance [START_REF] Amberkar | Diagnostic development for an electric power steering system[END_REF][START_REF] Peter | Electric power steering of the first step on the way to Steer by wire[END_REF].

The work described in this paper represents a part of development of new innovative brake handwheel concept, coupled to a Steer-by-wire system [START_REF] Messaoudene | Innovative brake handwheel concept for paraplegic[END_REF]. This device was intended initially for the paraplegic people, but it can be also used by the valid drivers considering these many advantages in safety and ergonomics. Two DC electrical engines are used in the designed prototype. First DC engine provides steering feedback on the handwheel by belt transmission. Braking action is obtained by exerting a pressure force on the handwheel equipped with a controllable braking feedback provided by second DC electrical engine with an active system by rack-pinion transmission. This solution allowed us to integrate dynamic brake assistance by using a command law which considers the displacement of the handwheel m braking action and its velocity.

The acceleration is simply obtained through an already available solution which consists of a circle behind the handwheel. The braking handwheel system developed in our laboratory is shown in Fig. 1. This concept allows the driver to brake and accelerate while keeping the hands on the handwheel in all driving situations. This also offers the advantage of a significant optimization in braking time and a more flexible and instinctive braking for a comfortable driving.

Fig. 1. Braking handwheel prototype

The use of X-by-wire technology in the automotive domain, offers many advantages in safety and driving assistance [START_REF] Inagaki | Yaw rate feedback braking force distribution control with control by-wire braking system[END_REF][START_REF] Bedner | A supervisory control to manage brakes and four-wheel-steer systems[END_REF][START_REF] Taneda | Design of actuator for active-rear-steer system[END_REF]. In the developed prototype, this technology permits a greater simplicity and precision of braking and acceleration actions by the driver during steering operations.

In this paper, we will be interested only on the modeling, simulation and experimental test of a Steer-by-wire system in the case of stopped vehicle.

In the first, we will describe the proposed dynamic model for Steer-by-wire system with hydraulic assistance. We will establish dynamic equations and commands laws for two DC electric engines included in the Steer-by-wire architecture. Simulation results will be given and commented in the end of this paper, by illustrating the handwheel feedback relevance obtained on the prototype used.

TIRE BEHAVIOR DURING STEERING IN STOPPED VEHICLE CASE

Tn the case of stopped vehicle, the front tires are sheared only by the resulting force coming from the driver torque on the handwheel and by the friction forces created in the pneumatic contact. When driver tries to steer the front wheels in stopped vehicle case, a road reaction force is applied in the tires and prevents front wheels to steering, by creating a sheared surface of an angle8 s (see Fig. 2). The front wheels begin to steer only when the side force on the tires reaches its saturation zone. Thus, this force becomes insufficient for opposing to the driver force on the handwheel [START_REF] Duysinx | Introduction a la mecanique des pneumatiques -comportement lateral[END_REF].

VVheel steering angle e w sheared surface Contrary to the moving vehicle case where the aligning torque on the front wheels tends to bring back the handwheel on the zero neutral position. It is known that when the vehicle is stopped, the neutral position of the handwheel is not null. Thus, it is a friction torque which will be felt on the handwheel. The neutral position of the handwheel depends on the equilibrium position of the front wheels, which is generally different to zero when the vehicle is stopped. When driver decreases torque applied on the handwheel, this one turns over to its initial position if the reaction force on the front wheels did not reach the saturation limit. As soon as the saturation limit is exceeded, the handwheel turns over to a more important position from where the steering action [START_REF] Mohellebi | Conception et realisation de systemes de restitution de mouvement et de retour haptique pour un simulateur de conduite[END_REF].

In order to restore on the handwheel the dynamic behaviour described, we will modeling the friction forces F i r on the tireroad contact by using the following equation [START_REF] Yih | Steer by wire : Implications for vehicle handing and safety[END_REF]:

With B, the rotation velocity of the sheared surface on the front tires. F , the side force developed on the tire-road contact, due to the load F,,{ applied on the front pneumatic. As shown in the Fig. 3, Steer-by-wire system is composed of two blocks: The first block provides to the driver a feedback on the handwheel, the second one ensures the steering act on the front wheels. A first electrical motor of steering is connected to the rack-and-pinion by reducer via the assistance valve. A second electrical motor of restitution, equipped with reducer, is connected to the handwheel [START_REF] Amberkar | A control system methodology for steer by wire systems[END_REF].

The signification of all the parameters used in modelling is shown on Table I.

Table 1. Signification of all used parameters The dynamics of the restitution block on the handwheel is described by the following equations.

Jh -B h J'T -B'T Jre -Bre Jw -Bw
Assessment of the efforts on the handwheel and the restitution motor gives:

wIth:

D = B _ � (3) jJJjh = -B/ )h + Td -KbhDb. B .. _ _ . Kbh b h R JreBre -BreBre + Tre + Db re R re
The state-representation of the restitution block on the handwheel is thus given by these two equations:

(4)
The steering block dynamics on the front wheels is described by the following equations:

.. The steering ratio R ra l i a corresponds to the transmission ratio between handwheel and front wheels around its pivot axes. It is often approximated by d / R knowing that the pivot arm p angle a . and trail plV angle CPt are relatively small [START_REF] Coudon | New reference model for Steer-by-wire applications with embedded vehicle dynamics[END_REF].

Replacing the expression of torsion angle on the assistance valve Dvalve ' one obtains:

.. . Kb darmKb " Bw (7)
J,eBse + B,eBse + R 2 " Bs e = T w + R R cos " cos a I" e se p "+'t piv

The steering system used is combined with a hydraulic assistance; modelling thus is very inspired of the real conventional steering. The efforts Assessment at the level of front wheels steering around its pivot axis gives:

Tt ra n s and T ass Is are respectively the driver torque transmitted to the rack-and-pinion and the assistance torque provided by the hydraulic system. These torques are applied to the rack and-pinion and obtained by the following expressions:

(9)

!(S va/ve ) is the law of hydraulic assistance valve and S the section of the assistance actuator. The product !(S v al v e ) S represents thus the assistance force applied to the rack-and pinion by the assistance actuator. Several dynamic models of the law valve function are proposed in [START_REF] Claeys | Commande laterale des vehicules automobiles avec prise en compte de la dynamique des actionneurs[END_REF], and [START_REF] Canudas | Fun to drive by feedback[END_REF] for static model. Tn our case, we linearized the law valve, around an operating point is adopted:

(

) 10
Where Ka is the slope of the function{ (S va l l' e ) at the selected operating point.

WIth : K K = -----'""--+ K S R a p (12)
Concerning the efforts applied on the sheared surface around the pivot axis of the front wheel, the dynamic equation is given by: J)J", = -Brde, + T" -Fi r ( t m + t p ) cos ¢t cos a piv [START_REF] Claeys | Commande laterale des vehicules automobiles avec prise en compte de la dynamique des actionneurs[END_REF] Where F r r is the friction force created on the pneumatic To obtain the instruction torque of the steering motor, an anolgy with the conventional steering system with a steering column is used. In this case, the steering engine fills the absence of the steering column, this condition should be satisfied: Using a similar procedure, we obtain the reference torque of the restitution engine on the handwheel. This torque represents an image of the external efforts due to the friction forces (in case of stopped vehicle) and lateral vehicle dynamics (in case of moving vehicle), which goes back to the handwheel via the steering column on a conventional steering system [START_REF] Wilwert | Influence des fautes transitoires et des performances temps reel sur la sfirete des systemes X-by wire[END_REF]. The reference torque is given by: T a dd is the additional torque which allows the adjustment of the handwheel feeling and its inertia, also to control the steering assistance on the handwheel according to the vehicle velocity by using a control law. This additional torque is given by this expression: , --Angle 55 0. 3 -----�------: -----� ------: -----� ------: ------� ------: ----: ----- As we can note it, these two curves show clearly the process of back move of the tire towards sheared surface when the front wheel undergoes an effort which is slackened abruptly. By analysing these curves, we distinguish three zones:

I I I I I I I I I 0. 05 -----� ------: ----� ------: ------� -----�------� -----�------� ----- , , , , , , , , , I I
1-The first zone corresponds to the phase where a force is applied on the front wheel. In this zone the tire begins to move without involving the movement of sheared surface, which remains motionless, because force resulting from the driver is lower than frictions force at the pneumatic contact.

2-In the second zone, the effort transmitted to the front wheel by the driver becomes higher than frictions force on the pneumatic contact. Sheared surface is pulled by the front wheel movement. The two angles follow the same movement.

3-In the third zone, the force applied to the front wheel is abruptly stopped. Tn this case, the wheel moves back towards the position of the sheared surface pulled by the reaction force due to the deformation of the tire in the pneumatic contact, while sheared surface remains in the same position.

Driver torque in right-left steering profile

Driver torque in right-left steering profile is used in this second simulation, showed on Fig. 7. The objective is to show the process of the steering hysteresis form. --------�-----;----------i ----------;-------- The driver torque according to the handwheel angle is shown in Fig. 8. We achieve a hysteresis diagram which characterizes perfectly real effect of feedback felt on the handwheel on real steering system when vehicle is stopped. -------'----------�-------------------: --------;-------- It's very important to notice that this proposed model was tested experimentally by using the developed prototype. With the adjustment of the various parameters we could obtain a handwheel feedback very near to that felt on the handwheel of a real steering system in the case of stopped vehicle. The effect of the handwheel back move to the position of the sheared surface is relevant.

6-

CONCLUSION

The complete replacement of the steering linkages by the x by-wire system poses problems related to the cost, standardization, power and reliability. The challenge of this technology lies in its capacity to detect and manage the breakdowns, by ensuring more safety of users. The constructors and researchers continue their developments in this field with an aim of increasing reliability and safety.

In this paper, we have first presented the context of the conducted research. Then, we have described the behaviour of the tire during the steering operations in the case of stopped vehicle. The process of the creation of a sheared surface at the wheel contact zone was explained, by showing the developed friction and reaction forces due to the deformation of the tire. The relevant part of this paper has been dedicated to the description of a proposed Steer-by-wire model in the case of stopped vehicle. Detailed dynamics study of this steering type with hydraulic assistance has been carried out, by giving all the relating dynamic equations.

A general Steer-by-wire architecture diagram has been explored. Simulation results have been presented. They have shown the coherence of the dynamic behaviour of the proposed model with hydraulic assistance in the case of stopped vehicle. In addition, the control law for the feedback torque on the handwheel has been tested using the developed prototype. The obtained results are very promising and the haptic feedback perceived on the handwheel is very realistic.

Fig. 2 .

 2 Fig. 2. Shearing of the front tire in steering for stopped vehicle

  the following expression, with f-J, f the adhesion coefficient of the front tires:3. USED STEER-BY-WTRE ARCHITECTURE & EFFORTS INTERACTIONSteer-by-wire architecture used in the modelling preserves the hydraulic assistance during steering actions[START_REF] Coudon | Algorithmes de commande pour Ie pilotage d'une direction decouplee[END_REF].

Fig. 3 .

 3 Fig. 3. Steering system with hydraulic assistance

  e h-e" e re -e se T d -Tss Tre -f'se R -R re I"e Inertia of handwheel -Its viscous friction Inertia of steering engine -Its viscous friction Inertia of handwheel restitution engine -Its viscous friction Inertia of front wheels -Its viscous friction Handwheel angle -Front wheels angle Angle of hand wheel restitution engine -Angle of steering engine Driver torque -Torque due to sheared surface Torque of hand wheel restitution engine -Torque of steering engine Reducer ratio of the restitution engine -reducer ratio of the steering engine R p -darm Pinion radius of the rack-and-pinion -Length of steering arm ¢t -a pi l' Trail angle of the wheels-Inclination of the pivot axis on the front wheels Kbh -Kbr Stiffness of the torsion bar on handwheel -Stiffness of assistance valve S " -R ra/ l n Torsion angle of bar -Global steering ratio between handwheel and wheels angle S -S va l ve Surface of assistance actuator -Torsion angle of the assistance valve

  . K J B = -B B + T ------'""--

  The torque T ., S is created by the reaction force ( e l' -e ,) c I in the sheared surface of the tire with an arm of distance ( t + t ) cos A.. cos a . . This torque takes the form: m p 'f/ t plV T ss = ( e lf -e s)CI ( t m + t p )C OS¢t cosa p l v (II) e., is the angle of sheared surface of the front tires. C I is the lateral stiffness of the front tires. t and t are respectively HI p mechanical and pneumatic trail of the tire. ¢t and a I' l l' are respectively the trail angle and pivot angle of the front wheels. The equation describing the steering dynamics of the front wheels becomes: e [ C. d J e [ darmK 1

  contact of the tire, shown previously. Brd is the viscous friction caused by the pneumatic contact on the road. The distance ( t + t ) cos do cos a . represents the arm of the HI p Y't (HF friction force applied in the sheared surface. The dynamic equation of the movement of sheared surfaces around pivot axis becomes: (19) J add and B a dd are parameters which allow the adaptation of handwheel felt. The suggested model for the Steer-by-wire dynamics in the Jwe, + Bri), + e , [ C 1 ( t m + t p )c os¢, cosa p ; v ] = (C Ar -F ! r )[( t m + t p )c os¢, cosa pi l' ] (14) case of a stopped vehicle is summarized in the form of a general diagram provided in the Fig. 4.

8 ":

 8 = � , e =? 8.w = 8" R.w se (15) By replacing this condition in the dynamic equation of the steering motor given previously, we obtain the formula of the torque command of steering engine according to the handwheel angle and the reference angle of front wheels e ,: e ! : .. . K T., c = J,c R , e e h + B,c R , e e h + R R p cos¢,cosa p ; v W (16) The reference angle of front wheels e ,: e 1 is obtained by replacing e , e = e h R ,c in the Equation (12).

!

  ivalve K bv + T a dd T re = -R rc = By replacing: e ,c = e h R , e ' one obtains:

Fig. 4 .

 4 Fig. 4. General diagram of the proposed steer-by-wire model Thereafter, some simulation results of the dynamic model suggested for a steering-by-wire system are presented.

4. 1 Fig. 5 .

 15 Fig. 5. Driver torque profile in teeth of saw form. The curves of the steering angle of front wheels and angle of sheared surface, obtained in simulation, are shown on Fig.6. 0. 3 5 ,-------,----:-::c:--,,--,-----,----, ---,---,------,----,

--

  Angle Wheel : ______ � -----� ------: ------:------: -----� ___ _

  1 5 -----i----+----;----+----i----i ------i-----+----+----J:: : ::::

Fig. 6 .

 6 Fig. 6. Steering angle and sheared surface angle

  ------�----------�--------� ----------: ----------� --------- ,, ,

- 4 - 6 Fig. 7 .

 467 Fig. 7. Driver torque in right-left steering profile

� 2 -t- 6 - 8 Fig. 8 .:

 2688 Fig. 8. Hysteresis form obtained in simulation. The centre shift of the obtained hysteresis curve Fig.8, is due to the reaction force created by the deformation of the tire in the contact zone, after slackening the handwheel at the end of each steering (right-left). The front wheels steering angle and the angle of sheared surface, are shown in the Fig. 9.

Fig. 9 .

 9 Fig. 9. Wheels and sheared surface angles We can notice that the obtained results in simulation for the dynamic Steer-by-wire model suggested in the case of a stopped vehicle are very conclusive and very coherent with the behaviour of a real steering system of vehicle with hydraulic assistance.


Member, IEEE Fakhreddine Ababsa
email: absbsa@iup.univ-evry.fr

Robust Camera Egomotion Estimation from 3D Straight Line-Based Environment model

Keywords: augmented reality, markerless tracking, extended kalman filter, Ransac. I

In this paper we present a new robust camera pose estimation approach based on 3D lines features. The proposed method is well adapted for mobile augmented reality applications We used an Extended Kalman Filter (EKF) to incrementally update the camera pose in real-time. The principal contributions of our method include first, the expansion of the RANSAC scheme in order to achieve a robust matching algorithm that associates 2D edges from the image with the 3D line segments from the input model. And second, a new powerful framework for camera pose estimation using only 2D-3D straight-lines within an EKF. Experimental results on real image sequences are presented to evaluate the performances and the feasibility of the proposed approach in indoor and outdoor environments

INTRODUCTION

Camera ego motion estimation is one of most challenging problem in computer vision. Several approaches based on natural features (corner points, planes, edges, silhouettes, etc.) in the scene have been developed last years. The mean idea of these techniques is to find correspondences between 2D features extracted from the image and 3D features defined in the world frame. The problem is then solved using 2D-3D registration techniques. Numerical nonlinear optimization methods like the Newton-Raphson or Levenberg-Marquardt algorithm are generally used for the minimization. Wuest et al. [START_REF] Wuest | Adaptive Line Tracking with Multiple Hypotheses for Augmented Reality[END_REF] present a model-based line tracking approach that can handle partial occlusion and illumination changes. The camera pose is computed by minimizing the distances between the projection of the model lines and the most likely matches found in the image. Drummond and Cipolla [START_REF] Drummond | Real-Time Visual Tracking of Complex Structures[END_REF] propose a novel framework for 3D model-based tracking. Objects are tracked by comparing projected model edges to edges detected in the current image. Their tracking system predicts the edge locations in order to rapidly perform the edge search. They have used a Lie group formalism in order to transform the motion problem into simple geometrics terms. Thus, tracking becomes a simple optimization problem solved by means of iterative reweighed least squares. Yoon et al. [21] present a model-based object tracking to compute the camera 3D pose. Their algorithm uses an Extended Kalman Filter (EKF) to provide an incremental pose-update scheme in a predictionverification framework. In order to enhance the accuracy and the robustness of the tracking against occlusion, they take into account the measurement uncertainties associated with the location of the extracted image straight-lines. Recently, Comport et al. [START_REF] Comport | Realtime markerless tracking for augmented reality: the virtual visual servoing framework[END_REF] propose a real-time 3D model-based tracking algorithm. They have used a visual servoing approach to formulate the pose estimation problem. A local moving edges tracker based on tracking of points normal to the object contours is implemented. In order to make their algorithm robust, they have integrated a M-estimator into the visual control law. Other approaches have also been applied where different features have been combined to compute the camera pose. Ababsa and Mallem [START_REF] Ababsa | Robust camera pose estimation combining 2D/3D points and lines tracking[END_REF] propose to combine point and line features in order to handle partial occlusion. They integrated a M-estimator into the optimization process to increase the robustness against outliers. Koch and Teller [START_REF] Koch | Wide-Area Egomotion Estimation from Known 3D Structure[END_REF] describe an egomotion estimation algorithm that takes as input a coarse 3D model of an environment . Their system uses a prior visibility analysis to speed initialization and accelerate image/model matching. Other approaches use Simultaneous Localization And Mapping (SLAM) to track the camera pose while building a 3D map of the unknown scene [START_REF] Davison | MonoSLAM: Real-Time Single Camera SLAM[END_REF] [START_REF] Williams | Automatic Relocalisation for a Single Camera Simultaneous Localisation and Mapping System[END_REF]. The main problem with most existing monocular SLAM techniques is a lack of robustness when rapid camera motions, occlusion and motion blur occur.

In this paper we present an original robust camera pose tracking using only straight lines and which differs from existing work. We propose to combine an EKF with a RANSAC scheme in order to achieve a robust 2D-3D lines matching. This gives an efficient solution for outliers rejection.

To our knowledge such solution has not been explored before. Furthermore, we have combined the 2D-3D lines correspondence constraints for object pose estimation, developed by Phong et al. [START_REF] Phong | Object Pose from 2D to 3D Point and Line Correspondences[END_REF], with an EKF in order to update recursively the camera pose. We have compared our results with classical approaches where pose estimation is solved using least square approaches [START_REF] Lowe | Fitting Parameterized Three-Dimensional Models to Images[END_REF][11] [START_REF] Lu | Fast and globally convergent pose estimation from video images[END_REF]. Our method requires no training phase, no artificial landmarks, and uses only one camera.

The rest of the paper is structured as fellows: In section II, we describe the camera pose estimation problem formulation when using straight lines, and we also give a complete implementation of the Extended Kalman Filter to update the camera pose recursively over the time using 2D and 3D lines features. In section III, we explain how we have expended the RANSAC scheme [START_REF] Fischler | Random Sample Consensus: A Paradigm for Model Fitting with Application to Image Analysis and Automated Cartography[END_REF] in order to achieve robust 2D-3D lines matching. In section IV, we show experimental results and evaluations, we discuss also the merits and the limitations of the proposed approach. Conclusion and future work are presented in section V.

II. CAMERA POSE ESTIMATION ALGORITH

In any Kalman Filter implementation, the system state is stored as a vector. In our algorithm the state is represented by the position and the orientation of the camera with respect to the world coordinate system. For computational we use a unit quaternion to represent the rotation. Thus, the state vector is given by:

[ ] z y x z y x t t t q q q q X 0 = (1)
where ( )

1 2 2 2 2 0 = + + + z y x q q q q .
We denote the camera state at time t by the vector t X . The EKF is used to maintain an estimate of the camera state X in the form of a probability distribution ( )

t t t Z X X P , 1 -
, where Zt is the measurement vector at time t. The Kalman filter models the probability distribution as Gaussian, allowing it to be represented by a covariance matrix Σ. In an extended Kalman filter, the non linear measurements and motion models are linearised about the current state estimate as Jacobian matrices.

Our algorithm follows the usual predict-refine cycle, whereby the state X is predicted at timestep t, and measurements are used to refine the prediction.

A. Time update

The time update model is employed in order to predict the camera pose at the following time step. In our case, the time update is simple because of the fact that we estimate the camera pose at each frame of the images sequence. Therefore the 3D camera pose between two successive frames changes very little. The time update equation is then given by :

1 - - ⋅ = t t X A X ( 2
)
Where A is 7×7 identity matrix.

The time update step also produces estimates of the error covariance matrix Σ from the previous time step to the current time step t. To perform this prediction we use the general update equation of the Kalman filter:

1 - - + ′ ⋅ Σ ⋅ = Σ t t t Q A A ( 3
)
Where Qt represents the covariance matrix of the process noise.

Σ reflects the variance of the state distribution.

B. Measurement model and estimate update

The measurement update model relates the state vector to the measurement vector. Since our goal is to estimate the camera pose using only straight lines, we will first describe the constraint equation which relates the state vector to the 3D model lines and their corresponding 2D image edges. We choose to base our technique on line features, rather than points, because this approach is relatively unexplored in the vision literature. We consider a pin-hole camera model and we assume that the intrinsic camera parameters are known. The world coordinate frame is a reference frame. All the 3D model lines are defined with respect to it. Let Li be a 3D line. Li is represented with the Cartesian coordinates of its two end-points i P 1 and i P 2 (see figure 1). The points i P 1 and i P 2 in world coordinates can be expressed in the camera frame as well :

+ ⋅ = + ⋅ = T i W P R i C P T i W P R i C P / 2 / 2 / 1 / 1 ( 4
)
Where the 3×3 rotation matrix R and the translation vector T describe the rigid body transformation from the world coordinate system to the camera coordinate system and are precisely the components of the camera state vector. We can see that the points

i C P / 1 , i C P / 2
and the center of projection OC are coplanar. i N is the unit vector normal to this plane. i N can be expressed in the camera coordinates frame as follows :

i C C i C C i C C i C C i P O P O P O P O N / 1 / 1 / 1 / 1 × × = (5)
Furthermore, a measurement input of the normal vector i N can be obtained from the image data. Indeed, image line matched with model line belongs also to the projection plane defined above. Let li be a 2D image line corresponding to the 3D line L i . In similar manner l i is represented by its two extremities

[ ] T i i i v u m 2 1 1 = and [ ] T i i i v u m 2 2 2 =
defined in the 2D image coordinates frame. The points i m 1 and i m 2 can be expressed in the camera coordinate frame as follows:

[ ] [ ] ⋅ = ⋅ = - - T i i i C T i i i C v u K m v u K m 1 1 2 2 1 / 2 1 1 1 / 1 (6)
Where the matrix K contains camera calibration parameters, such as focal length, aspect ration and principal point coordinates.

A measurement i n of the unit vector i N normal to the projection plane is thus given by (see figure 1):

i i i n n n 2 1 × = ( 7
)
Where

i C C i C C i i C C i C C i m O m O n et m O m O n / 2 / 2 2 / 1 / 1 1 = =
Combining equations ( 5) and ( 7), a measurement equation can be written, for each matching event L i l i :

( ) t t t v X h z + = ( 8 ) Where [ ] ( ) [ ]
= = = = T i z i y i x i t T i z i y i x i t N N N N X h n n n n z ( 9
)
v t represent the noise term in the measurement input with covariance R t . The noise is due to the uncertainty in the measured image position of the end points of the extracted 2D lines. The non linear function h(X) in measurement equation [START_REF] Williams | Automatic Relocalisation for a Single Camera Simultaneous Localisation and Mapping System[END_REF] relates the state to the measurement input. Three 2D-3D line correspondences are sufficient in theory to recover 6-DOF camera pose [START_REF] Dhome | Determination of the attitude of 3D objects from a single perspective view[END_REF] through in practice mores line may be required to increase accuracy.

The state estimate and covariance are refined after each feature measurement z t using the standard equation of the EKF as follows:

( ) ( ) ( ) - - - - - - - Σ ⋅ ⋅ - Σ = Σ - ⋅ + = ⋅ Σ ⋅ ⋅ ⋅ Σ = t t t t t t t t t t t T t t t T t t t H K X h z K X X R H H H K 1 (10)
Where H t is the Jacobian matrix defined by:

() - = ∂ ∂ = t X X t X X h H ( 11
)
The measurement update model is executed once a set of 2D-3D matched lines become available.

C. Iterated EKF

The standard EKF method does not consider errors due to the linearization of the non linear function h(X) in the vicinity of - t X . However, theses errors can lead to wrong estimates and/or divergence of the camera pose. Since the nonlinearity is only in measurement equation, the Iterated Extended Kalman Filter (IEKF) is the best technique to deal with it. The IEKF uses the same prediction equation as EKF, namely (2) and (3). The measurement update relations are replaced setting

- = t t X X 0
and doing iteration on :

( ) ( )

k t t k t k t k t t T k t t k t T k t t k t X X k t X h z K X X R H H H K X X h H k t - ⋅ + = + ⋅ Σ ⋅ ⋅ ⋅ Σ = ∂ ∂ = + - - - = 1 1 ) ( (12) For iteration number 1 , , 1 , 0 - = N k
. At the end of all iterations,

N t t X X =
. The covariance matrix is then updated based on N t X according to :

- - Σ ⋅ ⋅ - Σ = Σ t N t N t t t H K ( 13
)
The iteration could be stopped when consecutive values k t X and

1 + k t X
differ by less than a defined threshold III. ROBUST 2D-3D LINES MATCHING ALGORITH In this section we explain the expansion of the RANSAC scheme that we have developed in order to achieve a robust matching algorithm that associates 2D edges from the image with the 3D line segments from the input model, and without using any verification algorithm. Let 1. Randomly sample subsets of four {li ↔ Lj} pairs of 2D and 3D lines. In theory a minimum three pairs are of lines are sufficient to compute an accurate rigid transformation.

2. For each sample, compute the camera pose Π(R,T) using the IEKF algorithm described in section II.

3. Each candidate Π is tested against all the correspondences j i L l → by computing, in the camera frame, the angle between the normal vector i n (see figure 1) associated with the image line i l and the transformed line j L R ⋅ . If this match is wrong with respect the pose Π, then the co sinus of the angle should be significantly larger than zero. 4. We choose the pose Π which has the highest number of inliers, i.e the Π for which all the pairs are within a fixed angle threshold.

Hence, the obtained camera pose for the current image is robustly updated using only inliers of correspondences.

IV. EXPERIMENTAL RESULTS

The proposed camera pose estimation algorithm have been tested in real scenes and the registration accuracy was analyzed. The results are presented for both indoor and outdoor images (Figure 2). The indoor scene consists of a camera moving in an office room whereas the outdoor one corresponds to a moving camera pointing towards one frontage of a building. The frame rate of the recorded image sequences is about 25 frames/s and the resolution of the video images is 320×240 pixels. The 3D models of the office room and the building frontage are known, they are composed, respectively, of 19 lines and 120 line defined by the 3D coordinates of their end points within the world coordinates frame (Figure 3). In order to estimate the camera pose accuracy we defined, in the camera space, the registration error ξ. Given a set of correspondences between image edges and model segments the error ξ corresponds to the normalized square sum of the sinus of the angular disparities i α for each correspondence between image edge and the re projected model segments (Figure 1):

() = = × = = M i i i M i i N n M M 1 2 1 2 1 sin 1 α ξ (14)
Where M is the number of correspondences and α i the angle between the two planes spanned by the camera center, the observed image edge l i and the model segment i (see figure 1).

In our experiment we took M=10 correspondences. We have first considered that the data set has no outliers (100 % inliers or good matching) and we computed the registration error for several frames of the image sequence. Our algorithm was run on both the indoor and outdoor scenes. Similar results were obtained in both situations, Thus, for the indoor scenario, the mean error is about In the second experiment, we have evaluated the capacity of our robust algorithm to reject outliers in observed data. In our case, an outlier corresponds to a wrong feature matching between a 2D and a 3D line. For that, we have contaminated the input data set with different percentage of outliers, for both indoor and outdoor scenes, and have computed the corresponding registration error. The obtained results are summarized in table 1. As can be seen, our robust algorithm succeeded in all the cases to detect and delete the outliers. The camera pose is then estimated using the final data consensus which contains only the good 2D-3D lines correspondences. For example, in the worst case when 60% of the input data for the indoor scene (i.e. 6 lines correspondences among 10) are labeled as outliers, our algorithm was been able to identify the four inliers in the data. The camera pose returned using this inliers gives a registration error about 8.18×10 -5 . This result demonstrates the robustness and the accuracy of the proposed approach. Furthermore, we note that the number of trials needed to get the best solution increase with number of outliers (for example 240 trials for 60% of outliers for the indoor scenario and 56 trails for 50% of outliers for the outdoor scene). This means more processing time and will decrease the real time performance of the algorithm. 40% of outliers (15 trials on average) represents a good compromise. Figure 5 shows the projection of the models using the pose estimated by our algorithm for different frames of the indoor and outdoor images sequence and when run with different percentage of outliers. We can see that all lines are well aligned, and in this cases, the outliers have not affected the registration performance of our robust algorithm. Another advantage of our approach is its robustness to severe lines occlusion. Indeed, as the line constraint equation (see section II-B) for the camera pose parameters was developed in the case of "infinite image line". Any image points on the 2D line can be used to construct the corresponding projection plane. So, when partial occlusion occurs, it is enough to detect only small parts of the image edges to estimate the camera pose. In figure 6 we can see that several image edges are partially occluded (table and door for the indoor scene and because of the trees and the cars for the outdoor scene), in spite of that, the camera pose was successfully estimated.

We analyzed the processing time needed for camera pose estimation on a Pentium IV with 3GHz. All computations were performed in Matlab. The pose estimation process using IEKF does not take much time. Indeed, we have tuned the parameters of the IEKF in such manner so that it converges in few iterations (20 at the maximum). The processing time strongly depends only on the number of outliers in the current field of view. For example, the average time is about 28 millisecond per frame when having 40% of outliers in 10 input data for indoor scene. 3 milliseconds are used to estimate the camera pose with the IEKF and 25 milliseconds are measured for the time needed to reject outliers.

We have also evaluated the 3D camera localization accuracy. For that, we considered the outdoor scene, and we compared the computed camera trajectory obtained with our algorithm with the ground truth data. The ground truth for position are provided by a Trimble Pathfinder ProXT GPS receiver giving a submeter accuracy, The covered distance is about 30 meters. The ground truth for rotation was obtained using an Xsens MTx gyro mounted rigidly with the camera [32]. The results of the experiments are given in figure 6 which shows the error between the ground truth data and the recovered position and orientation of the camera for different frames of the image sequence. Table 2 reports the mean values and the standard deviation of the position and the orientation errors. It can be seen that our algorithm allows a good 3D localization performance, especially for orientation components.

V. CONCLUSION

In this paper, we proposed a new approach for 6-DOF camera localization based on matching between 2D image edges and 3D model segments. We performed a generic camera pose estimation framework based only on lines features using an Iterated Extended Kalman Filter. We also achieved significant improvements on robust 2D/3D lines matching scheme by adapting the well-know RANSAC algorithm to our application. The experimental results confirm the robustness of our approach against severe occlusion and outliers for both indoor and outdoor applications. We also showed the good performance of our algorithm to localize a moving camera in 3D environment. Future work will be devoted to extend our approach by using other outdoor sensors (e.g. an inertial sensor and a GPS). Thus, the system could be used for navigation and localization in large-scale outdoor environments. An hybrid algorithm will the fuse the data provided by the three sensors in order to refine the camera pose.

Figure 1 .

 1 Figure 1. Projection plane. The model line, its projection onto the image and the center of projection OC are coplanar.

Figure. 2 .Figure. 3 .

 23 Two frames from the recorded images sequence. (a) indoor scene. (b) outdoor scene. 3D models of indoor and outdoor scenes used for experiments. (a) the office room model. (b) The frontage building model.

  shows the projection of the office model using the camera pose estimated by our algorithm, and as can see it is quite skewed. All the lines are fairly well aligned.

  . Projection of the 3D models using the final camera pose estimate when the input data has no outliers. (a) indoor scene (b) outdoor scene

  Outliers = 37% (d) outliers = 50% Figure. 5. Camera pose estimation results

Figure. 6 .

 6 Figure. 6. Position and orientation accuracy with respect to ground truth

  {}

	i , = image and { } l i 1 ,..., j L j , =	N ,..., be a set of 2D edges extracted from the M 1 a set of 3D model lines. Our robust
	lines matching algorithm is summarized as follows :

TABLE 2

 2

			3D LOCALIZATION PERFORMANCE
		Position errors (m)		Rotation errors (deg)
		ΔX	ΔY	θx	θy	θz
	Mean	2.26	1.42	0.81	2.61	1.38
	Std	0.88	0.77	0.66	0.41	0.75


Gilbert Bergna
email: gilbert.bergnadiaz@supelec.fr

Member, IEEE Matteo Pirro

Erik Berne
email: erik.berne@edf.fr

Amir Arzandé

Marta Molinas Cabrera

Member, IEEE Romeo Ortega

IEEE Fellow'99

Modular Multilevel Converter Passivity-Based PI Control suited for Balanced and Unbalanced Grid Conditions

MMC Passivity-based PI Controller: A novel nonlinear global tracking controller was proposed by extending the regulation case. It consists in a linear PI that regulates two constructed signals with respect to which the incremental model of the MMC becomes passive. Estimation of x ⋆ : A closedloop estimator was proposed to estimate the MMC variables in steady state needed as input for the non-linear controller. Such estimator consists on a virtual MMC energy model and both of the differential current reference equations that yielded from mathematical optimization using Lagrange multipliers.

I. INTRODUCTION

Voltage Source Converter (VSC)-based High Voltage Direct Current (HVDC) links [START_REF] Flourentzou | Vsc-based hvdc power transmission systems: An overview[END_REF], are currently the key technological solution for bulk power transmission over long distances. Moreover, multi-terminal HVDC connections have the potential of shifting Europe's power generation, from conventional to renewable energy sources. Such SuperGrid will serve as a transcontinental highway for renewable energy [START_REF] Ahmed | Prospects and challenges of future hvdc supergrids with modular multilevel converters[END_REF], allowing the geographical smoothing effects minimize the disadvantages inherent to the intermittent nature of renewable generation. The Modular Multilevel Converter introduced in [START_REF] Lesnicar | An innovative modular multilevel converter topology suitable for a wide power range[END_REF] will be the power electronic converter suited for such task [START_REF] Ahmed | Prospects and challenges of future hvdc supergrids with modular multilevel converters[END_REF], presenting several advantages with respect to its predecessors such as its high modularity, scalability and lower losses, although bringing challenging control and stabilization issues that must be addressed.

There has been a clear tendency of controlling the modular multilevel converter (MMC) in its natural "abc" reference frame, instead of the commonly used "dqo" or "αβo" frame that results from applying the well known Park or Clarke transformations. This interest arises in order to be able to regulate in a simple and direct way the energy stored in each phase of the MMC, unlike the three-level VSC that presents only solid state switches unable to store energy. The control strategy proposed in [START_REF] Bergna | A generalized power control approach in abc frame for modular multilevel converter hvdc links based on mathematical optimization[END_REF] has proven to be robust thanks to its closed-loop nature and its independence to the system parameters, as well as it has presented a good dynamic performance. Yet, such strategy is based on linear controllers, and since the system is highly non-linear, makes it impossible to ensure its stability.

Global asymptotic stability issues have been investigated by means of Lyapunov's proof in [START_REF] Antonopoulos | Global asymptotic stability of modular multilevel converters[END_REF], [START_REF] Harnefors | Global asymptotic stability of modular multilevel converters with measurement lag and circulating-current control[END_REF]. Indeed, asymptotic stability of the MMC under open loop and closed loop control was proven, yet the validity of such proof only holds for the system under study. This means that even if global asymptotic stability was proven for two different MMC converters, the stability of the interconnection between these two systems is not ensured, and to prove it using Lyapunov's stability criteria one must consider the new and more complex system as a whole. Clearly, since we are aiming at a multi-terminal HVDC systems which are supposed to grow "organically" in Europe, the final state of such SuperGrid is unknown since it will be continuously expanding. With the complete picture in mind, it seems almost impossible, or at least of a great degree of complexity to ensure global asymptotic stability of the ever growing multi-source SuperGrid using a global approach. Therefore a modular control approach for ensuring stability of the global interconnected system is needed. Last but not least, the required technique must be simple to implement regardless of the complexity of the proof.

Passivity is a characterization of the system behavior based on energy, and it can be used to asses the stability of a single system [START_REF] Ortega | Passivity-based Control of Euler-Lagrange Systems: Mechanical, Electrical and Electromechanical Applications[END_REF]. Most importantly, the negative feedback combination of two passive systems results in a passive and stable interconnection [START_REF] Mccourt | Stability of networked passive switched systems[END_REF], [START_REF] Stramigioli | Passive output feedback and port interconnection[END_REF]. Therefore, the idea of making complex multi-terminal MMC-HVDC systems based on local passivity controllers as a modular approach for ensuring global asymptotic stability of the global and much larger system as a whole may be considered as the philosophical contribution of this article.

Moreover, it is possible to use a PI controller to render the system globally asymptotically stable via passivity, simplifying the implementation as was done in [START_REF] Jayawardhana | Passivity of nonlinear incremental systems: Application to pi stabilization of nonlinear rlc circuits[END_REF]- [START_REF] Hernandez-Gomez | Adaptive pi stabilization of switched power converters[END_REF]. Such controllers are known to be simple, robust and widely accepted by practitioners. However, the aforementioned work requires that the desired state variables in steady state that become the references values (x ⋆ ) of the system should be constant values, in order to solve this regulation problem. As mentioned earlier, the aim is to control the MMC in the "abc" frame which implies that the reference state variables of the system, that is to say, the differential and grid currents (i di f f ⋆k and i v⋆k ), and the upper and lower capacitor arm voltages sum (u Σ cu⋆k and u Σ cl⋆k ) are not constants, turning the regulation problem of [10]- [START_REF] Hernandez-Gomez | Adaptive pi stabilization of switched power converters[END_REF] into a tracking problem instead.

Hence the extension of the passivity-based PI controller of [START_REF] Jayawardhana | Passivity of nonlinear incremental systems: Application to pi stabilization of nonlinear rlc circuits[END_REF]- [START_REF] Hernandez-Gomez | Adaptive pi stabilization of switched power converters[END_REF] from the regulation case to the tracking scenario is applied to the MMC as one of the main contributions of this work. However, the success of the implementation of this new control strategy strongly depends on the estimation of the steady state signals x ⋆ that are to be tracked. Such crucial issue is thoroughly discussed in this article, and a closed loop estimator using a virtual energy model of the MMC along with two different differential current control references that yielded from an analytical optimization problem based on Lagrange multipliers [START_REF] Bergna | A generalized power control approach in abc frame for modular multilevel converter hvdc links based on mathematical optimization[END_REF], [START_REF] Bergna | A generalized power control approach in abc frame for modular multilevel converters based on lagrange multipliers[END_REF] is proposed, capable of controlling the MMC in a phase-independent fashion or for preventing power oscillations caused in the unbalanced AC grid to flow to the DC side of the MMC.

The remainder of this article is organized as follows. The passivity-based PI controller is presented in section II. In order to successfully apply the control strategy, one must be able to estimate the steady state trajectories of the MMC variables. In IV is discussed how to calculate these tracking references for phase independent control and for constant DC power during unbalanced grid conditions. Simulations results are performed and their results are given in Section V. Finally, conclusions in Section VI complete the brief.

II. CONTROL THEORY

A. Global Tracking Problem

Consider the bilinear system

ẋ(t) = Ax(t) + d (t) + m ∑ i=1 u i (t)B i x(t) (1)
where 1 x, d ∈ R n are the state and the (measurable) disturbance vector, respectively, u ∈ R m , m ≤ n, is the control vector, and A, B i ∈ R n×n are real constant matrices.

Given an admissible, differentiable trajectory, that is a function

x ⋆ : R + → R n verifying ẋ⋆ =Ax ⋆ + d + m ∑ i=1 u ⋆ i B i x ⋆ (2)
for some control signal u ⋆ : R + → R m . Find, if possible, a dynamic state-feedback controller of the form

ż = F(x, x ⋆ ) (3) u = H(x, x ⋆ ), (4)
where F : R n × R n → R q , q ∈ Z + , and H : R n × R n → R m , such that all signals remain bounded and

lim t→∞ [x(t) -x ⋆ (t)] = 0, (5)
1 For brevity, in the sequel the time argument is omitted from all signals.

for all initial conditions (x(0), z(0

)) ∈ R n × R q .
In this paper we characterize a set of matrices {A, B i } for which it is possible to solve the aforementioned global tracking problem with a simple linear PI controller. The class is identified via the following LMI.

Assumption 1: ∃P ∈ R n×n , P = P ⊤ > 0 such that sym(PA) ≤ 0 (6) sym(PB i ) = 0, (7)
where sym : R n×n → R n×n computes the symmetric part of the matrix.

To simplify the notation in the sequel we define the positive semidefinite matrix

Q := -sym(PA). (8)

B. Passivity of the Bilinear Incremental Model

Lemma 1: Consider the system (1) verifying the LMI of Assumption 1 and an admissible trajectory x ⋆ . Define the incremental signals (

•) := (•) -(•) ⋆ ,
and the output function

y := C (x ⋆ )x (9)
where the map C : R n → R m×n is defined as

C :=    x ⊤ ⋆ B ⊤ 1 . . . x ⊤ ⋆ B ⊤ m    P.
The operator ũ → y is passive with storage function

V ( x) := 1 2 x⊤ P x. ( 10
)
Thus, the system verifies the dissipation inequality V ≤ ũ⊤ y.

C. A PI Global Tracking Controller

From Lemma 1 the next corollary follows immediately.

Corollary 1: Consider the system (1) verifying the LMI of Assumption 1 and an admissible trajectory x ⋆ in closed loop with the PI controller

ż = -y u = -K p y + K i z + u ⋆ (11)
with output [START_REF] Stramigioli | Passive output feedback and port interconnection[END_REF] and

K p = K ⊤ p > 0, K i = K ⊤ i > 0.
For all initial conditions (x(0), z(0)) ∈ R n ×R m the trajectories of the closedloop system are bounded and lim t→∞ y a (t) = 0, [START_REF] Castaños | Proportional plus integral control forset-point regulation ofaclass ofnonlinear RLC circuits[END_REF] where the augmented output y a : R + → R m+n is defined as

y a := C Q x. Rcharge Lcharge va vdc d (N) ±v adc(N) d (2) ±v adc(2) d (1) ±v adc(1) h (N) ±v ahc(N) h (2) ±v ahc(2) h (1) ±v ahc(1) R,L R,L ... ... Phase-a iLa iUa iva v ad(N) v ad(2) v ah(2) v ah(1) it v ki( j) Ck Sub-MODULE h (k)
C Q ≥ n, ( 13
)
then global tracking is achieved, i.e., (5) holds.

Remark 1: Notice that the matrix C depends on the reference trajectory. Therefore, the rank condition (13) identifies a class of trajectories for which global tracking is ensured.

Remark 2: Condition ( 13) is sufficient, but not necessary for state convergence. Indeed, as shown in [START_REF] Van Der Schaft | L 2 -Gain and passivity techniques in nonlinear control[END_REF], global tracking is guaranteed if y a is a detectable output for the closed-loop system. That is, if the following implication holds

y a (t) ≡ 0 =⇒ lim t→∞ [x(t) -x ⋆ (t)] = 0.

III. MMC APPLICATION

The model of the MMC converter (fig. 1) is

idi f f = - R L i di f f - 1 4L u CΣ u Σ - 1 4L u C∆ u ∆ + v dc 2L iv = - R ′ L ′ i v - 1 4L ′ u C∆ u Σ - 1 4L ′ u CΣ u ∆ - V pcc L ′ uCΣ = 1 C ′ i di f f u Σ + 1 2C ′ i v u ∆ uC∆ = 1 2C ′ i v u Σ + 1 C ′ i di f f u ∆ (14)
Four state variables are considered: 1) The differential or circulating current of the MMC i di f f , 2) the grid or load current i v , 3) the sum between the upper and lower capacitor voltages and4) the difference between them u C∆ = u CUu CL = u ∆ C . Then, using the notation of (1) we define

u CΣ = u CU + u CL = u Σ C ,
A =     -R L 0 0 0 0 -R ′ L ′ 0 0 0 0 0 0 0 0 0 0     , d =     v dc 2L - V pcc L ′ 0 0     B Σ =     0 0 -1 4L 0 0 0 0 -1 4L ′ 1 C ′ 0 0 0 0 1 2C ′ 0 0     , B ∆ =     0 0 0 -1 4L 0 0 -1 4L ′ 0 0 1 2C ′ 0 0 1 C ′ 0 0 0     ,
and the matrix P of assumption ( 1) is

P =     2L 0 0 0 0 L ′ 0 0 0 0 C ′ 2 0 0 0 0 C ′ 2     . It follows that, defining x ⊤ := [i di f f i v u CΣ u C∆ ]
⊤ , the passive output of the system is

y MC = x ⊤ ⋆ B ⊤ Σ x ⊤ ⋆ B ⊤ ∆ Px = 1 2 i di f f ⋆ u CΣ -i di f f u CΣ⋆ + 1 2 i v⋆ u C∆ -1 2 i v u C∆⋆ i di f f ⋆ u C∆ -i di f f u C∆⋆ + 1 2 i v⋆ u CΣ -1 2 i v u CΣ⋆ . (15)
Detectability is given by the rank of the matrix

C Q =        -u CΣ⋆ -1 2 u C∆⋆ i di f f ⋆ 1 2 u C∆⋆ -u C∆⋆ -1 2 u CΣ⋆ 1 2 i v⋆ i di f f ⋆ 2R 0 0 0 0 R ′ 0 0 0 0 0 0 0 0 0 0       
which satisfies condition (13) whenever i 2 di f f ⋆ -1 4 i v⋆ u C∆⋆ = 0 which is satisfied in practice.

IV. REFERENCES CALCULATION

A. Generation of references in an Open Loop fashion for Phase-Independent Control

The grid current reference equilibrium point is imposed based on the desired voltage reference of the load. For the system under study with no grid connection (i.e, V pcc = 0) such reference is expressed as follows:

i v⋆ = |e v⋆ | R ′2 + (ωL ′ ) sin ωt -tan -1 ωL ′ R ′ (16)
Where e v⋆ is the sinusoidal voltage that one wishes to apply to the load, with an amplitude that can vary from 0 to V dc 2 . In addition, R ′ = R 2 + R load and L ′ = L 2 + L load . The differential or circulating current reference i di f f ⋆ is estimated by assuming a small internal resistance R:

i di f f ⋆ ≈ |e v⋆ i v⋆ | V dc ( 17
)
Since i di f f ⋆ is constant (for this case), the voltage that drives it may be expressed simply by

u di f f ⋆ = Ri di f f ⋆ .
Also, it is possible to calculate the fluctuations of the sum w Σ⋆ and difference w ∆⋆ of the capacitive energy stored between the upper and lower arms of the MMC (see [START_REF] Harnefors | Dynamic analysis of modular multilevel converters[END_REF] for more details on such equations). Thus,

∆w Σ⋆ = t 0 -e v⋆ i v⋆ + V dc -2u di f f ⋆ i di f f ⋆ dt ∆w ∆⋆ = t 0 i v⋆ 2 V dc -2u di f f ⋆ -2e v⋆ i di f f ⋆ dt ( 18
)
Remark 3: A high-pass filter is needed for the term that is being integrated in equation [START_REF] Gnanarathna | Efficient modeling of modular multilevel hvdc converters (mmc) on electromagnetic transient simulation programs[END_REF] to leave out the power error caused by neglecting R in equation [START_REF] Angquist | Open-loop control of modular multilevel converters using estimation of stored energy[END_REF]. The average value of w Σ⋆ and w ∆⋆ is given as a reference from the user; typically, W Σφ⋆ = C N V 2 dc , and W ∆φ⋆ = 0. Using the energy estimation, it is now possible to calculate the upper and lower arm voltages.

u CU⋆ = N C W Σφ⋆ + ∆w Σ⋆ + W ∆φ⋆ + ∆w ∆⋆ u CL⋆ = N C W Σφ⋆ + ∆w Σ⋆ -W ∆φ⋆ + ∆w ∆⋆ (19)
The remaining state variables u Σ⋆ and u ∆⋆ may now be calculated as

u CΣ⋆ =u CU⋆ + u CL⋆ u C∆⋆ =u CU⋆ -u CL⋆ ( 20
)
In addition, the upper and lower insertion indexes are calculated by

n u⋆ = V dc 2 -e v⋆ -u di f f ⋆ u cu⋆ n l⋆ = V dc 2 + e v⋆ -u di f f ⋆ u cl⋆ ( 21
)
The control at the equilibrium point is therefore defined as

u Σ⋆ =n u⋆ + n l⋆ u ∆⋆ =n uk⋆ -n l⋆ (22)
This estimation technique has been implemented successfully along with the so-called open loop control strategy in [START_REF] Antonopoulos | Global asymptotic stability of modular multilevel converters[END_REF], [START_REF] Angquist | Open-loop control of modular multilevel converters using estimation of stored energy[END_REF]. Nonetheless, its effectiveness relies on the approximation made in [START_REF] Angquist | Open-loop control of modular multilevel converters using estimation of stored energy[END_REF] where the power dissipated in the internal resistance of the MMC is neglected. For high values of R, this assumption does not hold any more [START_REF] Harnefors | Global asymptotic stability of modular multilevel converters with measurement lag and circulating-current control[END_REF], and the estimating technique will no longer be useful. A solution for this case is proposed in the following sections.

B. Generation of references with a closed loop estimator

If the converter has a non-negligible internal resistance R > 0, a mismatch between the average values of powers e vk i vk and v dc i di f f k will take place, making approximation [START_REF] Angquist | Open-loop control of modular multilevel converters using estimation of stored energy[END_REF] no longer valid. To overcome this issue, a closed loop estimator is proposed to calculate x ⋆ and u ⋆ , and is depicted in Fig. 2. The estimation technique is based on the application of equation ( 23): Σk will "add" the missing power that was neglected by the previous estimation technique. The difference between the proposition of [START_REF] Bergna | A generalized power control approach in abc frame for modular multilevel converter hvdc links based on mathematical optimization[END_REF], [START_REF] Bergna | A generalized power control approach in abc frame for modular multilevel converters based on lagrange multipliers[END_REF] and what is being proposed here, is that equation (23) will be applied to a virtual MMC instead of the real converter; that is to say that it will not participate directly in the control of the converter, yet it will generate the references needed for the control by means of the MMC state equations. The procedure is described in the following lines.

i di f f ⋆ = P re f Σk + (1 -α) P vk v 2 dc v dc + αP vk v dc + -P re f ∆k
For a passive grid (V pcck = 0) the grid/load current is defined as in 16. If the system is conected to a three-phased active grid (V pcck = 0) then the grid current reference i v⋆ is defined by equation ( 24):

i v⋆ = P re f ac v + pcc 2 + kp • v - pcc 2 • v + pcck + kp • v - pcck ( 24
)
Equation ( 24) calculates the grid current in steady state (i v⋆ ), as a function of the power reference established by the secondary control P re f ac , and the positive and negative voltage measurements of the point of common coupling (v + pcc and v - pcc ) from the real MMC, and a constant kp defined by the user.

In order to calculate the differential current reference in steady state (i di f f ⋆ ) equation 23 is used along with a virtual energy model of the MMC. Such model contains the MMC energy sum and difference dynamics ( ẇΣ and ẇ∆ )

w Σ⋆k = t 0 -e v⋆k i v⋆k + V dc -2u di f f ⋆k i di f f ⋆k dt w ∆⋆k = t 0 i v⋆k 2 V dc -2u di f f ⋆k -2e v⋆ i di f f ⋆k dt (25) (26)
with e v⋆k and u di f f ⋆k calculated as

u di f f ⋆k =Ri di f f ⋆k + L d dt i di f f ⋆k e v⋆k =R ′ i v⋆k + L ′ d dt i v⋆k +V pcck (27)
The average values of the w Σ⋆ and w ∆⋆ are required to calculate P re f Σ⋆ and P re f ∆⋆ expressed in (28) that are in turn needed for the calculation of i di f f ⋆ by (23).

P re f Σ ⋆ k = kp Σ W re f Σk -w Σ k + ki Σ W re f Σk -w Σ k dt P re f ∆ ⋆ k = kp ∆ W re f ∆k -w ∆ k + ki ∆ W re f ∆k -w ∆ k dt (28)
In order to so whilst obtaining good dynamic performance, the average values of w Σ⋆ and w ∆⋆ are calculated using:

ŵΣ⋆ =w Σ⋆ -w Σ⋆α ŵ∆⋆ =w ∆⋆ -w ∆⋆α ( 29
)
Where w Σ⋆α and w ∆⋆α are band-pass-filtered signals produced by a second order generalized integrator (SOGI)quadrature signal generator that is acting as a Notch filter. Furthermore, u cu⋆ and u cl⋆ are calculated via (30)

u CU⋆ = N C [w Σ⋆ + w ∆⋆ ] u CL⋆ = N C [w Σ⋆ -w ∆⋆ ] (30)
Finally, u cΣ and u c∆ are calculated as in 20. This concludes the estimation of x ⋆ . To obtain u ⋆ ; i.e., the control in steady state, equations 21 and 22 are required once more.

C. Generation of references for constant DC power under unbalanced grid conditions

Nevertheless, i di f f ⋆ determined by equation ( 23) is not able to successfully cope with unbalanced grid conditions, since it is the result of an analytical mathematical optimization that does not consider the relationship between the phases of the MMC. To be able to successfully handle unbalances, by preventing power fluctuations from the unbalanced AC side of the converter to pass its DC side, i di f f ⋆ is calculated using equation (31): As was done for the previous case, equation (31) must be combined with the virtual energy model of the MMC expressed in (26). This is since (31) requires closing the energy regulation loop with three the PI controllers represented by P re f Σ⋆k , P re f ∆⋆ and, for the direct power assignment case, P dc . 1) The Controller: To summarize, the controller equations are:

i di f f ⋆k = P re f Σ⋆k + (1 -α) P v⋆k v 2 dc v dc + αP v⋆k v dc + -P re f ∆⋆k e v⋆k
żMC = -y MC u MC = -K p y MC + K i z MC + [u Σ⋆ u ∆⋆ ] ⊤ y MC = 1 2 i di f f ⋆ u CΣ -i di f f u CΣ⋆ + 1 2 i v⋆ u C∆ -1 2 i v u C∆⋆ i di f f ⋆ u C∆ -i di f f u C∆⋆ + 1 2 i v⋆ u CΣ -1 2 i v u CΣ⋆
With reference variables i v⋆ defined in (24), i di f f ⋆ in (23) for phase independent control and in (31) for constant DC power control; u CΣ⋆ and u C∆⋆ are defined by replacing (30) in [START_REF] Wang | Pliant active and reactive power control for grid-interactive converters under unbalanced voltage dips[END_REF]. In addition, u Σ⋆ and u ∆⋆ are derived in (21) and (22).

V. RESULTS

To validate the proposed control, simulations have been carried out.

i di f f k & i vk (Amps) u Σ c u,l (Volts) 0 5 • 10 -2 0.1 -4 -2 0

A. Phase Independent Control

A single-phase MMC simulation scenario has been set up in Matlab/Simulink using a high efficiency model [START_REF] Gnanarathna | Efficient modeling of modular multilevel hvdc converters (mmc) on electromagnetic transient simulation programs[END_REF] to test the validity of the control. The considerations are the following : the converter has 2N = 10 submodules, N in each arm (upper and lower). The input DC voltage is V dc = 140V , the reference voltage e v⋆ has an amplitude of V dc 2 and a frequency set to f = 50Hz. The frequency of the balancing algorithm [START_REF] Harnefors | Dynamic analysis of modular multilevel converters[END_REF] that balances the N capacitor voltages is set to f = 20kHz. The internal capacitance, resistance and inductance are respectively set to C = 3.3mF, R = 12Ω and L = 10mH. The load resistance and inductance values, respectively, are R = 12Ω and L = 40mH. Figure 3-a) shows the 2N voltage trends of the MMC capacitors, whereas in 3-b) is depicted the AC output voltage multilevel waveform of the converter. Furthermore, the state references x ⋆ as well as the measured state variable are depicted in Fig. 4. As can be seen from Fig. 4-a) and Fig. 4-b), the i di f f and i v converge to i di f f ⋆ and i v⋆ . Nonetheless in Fig. 4-c) a non-negligible error between u Σ c u,l ⋆ and u Σ c u,l of ≈ 1Volt can be observed. Fortunately, u Σ c u,l is the sum all the N capacitors in the upper or lower arm, thus the error of the individual capacitor voltage u c i ≈ 1/5Volts. Such error occurs since there is still a remaining power mismatch, even after reducing it using the proposed closed loop estimator along with the virtual energy model. This occurs since the high frequency current harmonics caused by the modulation, and eventually by the saturation of the insertion indexes, are not considered in such model. These current harmonics produces losses when passing through the MMC resistances, and cause the error in the estimation.

B. Constant DC Power Control

In this section, a three-phased MMC converter connected to an active grid has been simulated under unbalanced conditions in order to test the proposed control along with the closedloop estimator for constant DC power control. The simulation scenario was done by establishing the primary power reference of the system at the point of common coupling via P re f ac , as well as establishing it at the DC terminals of the MMC via P re f dc . The MMC parameters are given in table I. This series of simulation analyzes the performance of the controller when at t = 0.1s the voltage at the point of common coupling of phase "a" is reduced to 5% of its original value and remains at that value for a total time of 0.3s, before returning to its original level at t = 0.4. This is depicted in Fig. 5.

In Fig. 6-a) are depicted the power at the point of common coupling P ac in blue and the power at the DC terminals of the MMC P dc in red. It can be seen that during the fault P ac starts fluctuating at twice the grid frequency. Nonetheless, such oscillations are blocked by the energy storing elements of the MMC, and do not appear at the DC side of the converter. This holds independently of how the primary power reference of the system is being established. However, for the primary power reference is established at the DC terminals of the MMC, it seems that the power oscillations are most effectively decoupled as even transient perturbations have been removed. Nonetheless, this case imposes a more elevated energy and voltage stress to its capacitors as can be seen in Fig. 6-b) and 6-c), and have to be oversized to withstand this scenario.

In Fig. 6-c), 6-d) and 6-e) are shown the waveform trends of the references and the state of the system (x ⋆ and x): the differential currents, the grid currents, and the capacitor voltages, respectively. A familiar result comes to mind: the grid and circulating currents are regulated quite well, whereas a stationary regulation error is found between reference and measurements of the capacitor voltages. The error is approximately of the order of 2%. Here again, it is caused by the power mismatch error in the estimation. Furthermore, the grid currents stay balanced in both cases since kp = 0 is used in (24) [START_REF] Rodriguez | Flexible active power control of distributed power generation systems during grid faults[END_REF], [START_REF] Wang | Pliant active and reactive power control for grid-interactive converters under unbalanced voltage dips[END_REF].

VI. CONCLUSIONS

In this brief, a passivity-based PI controller has been applied to the modular multilevel converter in the ABC frame. Global asymptotic stability is ensured with a simple PI regulating the MMC passive outputs to zero. The obtained stability results are global and hold for all positive definite gains of the PI. The performance of the controller was tested by means of some realistic simulations for the MMC. The converter showed good performance despite the estimation error caused by neglecting the high order harmonics and physical limits including control saturation. The proposed estimation technique is capable of considering the power losses caused by the fundamental harmonics of each current in the internal resistance of the MMC, usually neglected in "open loop" approaches. Furthermore, it can be used to generate the necessary references x ⋆ necessary for preventing power oscillations caused by an unbalanced grid from appearing at the DC terminals of the MMC. However, the estimation is model-based, which is usually a drawback during the practical implementation due to the parameters uncertainty. Last but not least, passivity offers a local approach to ensure stability of complex MMC-based multi-terminal HVDC systems

Figure 1 .

 1 Figure 1. a) One-leg topology of a modular multilevel converter

Figure 2 .

 2 Figure 2. Closed loop estimator for phase independent control

  calculates the reference of the differential current of the MMC. Such equation contains two PI controllers represented by the terms P re f Σk and P re f ∆k that regulate the average values of the capacitive energy sum and difference stored in the arm of the converter ( w Σk and w ∆k ). In the case the MMC would have a non-negligible internal resistance R, the PI controller P re f

Figure 3 .

 3 Figure 3. MMC 2N Individual Capacitor Voltages, b) MMC output voltage

Figure 4 .

 4 Figure 4. Phase-independent control of a single phased MMC: a) Differential and Grid Current; b) MMC Sum of capacitor voltages by arm.

Figure 6 .

 6 Figure 6. Power and energy trends: a) Power at the dc terminals (in red) and at the PCC (in blue), b) MMC 3-phase capacitive energy, and c) MMC individual capacitor voltages d) differential currents and e) grid currents

Table I .

 I MMC PARAMETERS

			Parameter		Value
			DC Voltage	150 Volts
		Number of SM per arm	5
			SM Nominal Voltage	30 Volts
			SM Capacitance	3.3 mF
			Arm Inductance	20 mH
			Arm Resistance	6 Ω
			Load Inductance	20 mH
			Load Resistance	6 Ω
		50
	V pcc (Volts)	0
		-50
		0	0.1	0.2	0.3	0.4	0.5
				time (s)
	Figure 5. Voltage at the point of common coupling