Houari Bettahar

Olivier Lehmann

Cédric Clévy

Nadège Courjal

Philippe Lutz

Photo-Robotic Extrinsic Parameters Calibration of 6-DOF Robot for High Positioning Accuracy

Keywords: Nano-positioning, Accuracy, Photo-robotic, robot calibration, Photonics

Hybrid integration of individual photonic elements appears particularly promising. It indeed provides high performances, proposes new optical functionalities and products and exploits new propagation modes of light beams. This approach requires an accurate multi Degree-Of-Freedom (DOF) positioning of the individual photonic elements. Hence, the inaccurate multi-DOF measurements and robots control are the main locks to overcome, notably at the micro-scale. For this sake, an original photo-robotic approach has been proposed. It relies on multi-DOF robot's motion associated with the use of 1-D Fabry-Perot interferometry measurement to achieve multi-DOF pose measurement. This original photo-robotic approach notably integrates the issue of extrinsic parameters calibration through a proposed photo-robotic calibration model. Experimental investigations demonstrate that a translational positioning accuracy of 390 nm was obtained. This photo-robotic approach has especially been applied to achieve the 6-DOF positioning of an optical lamella relative to an optical fiber with high accuracy that also conducts to maximum optical reflectivity.

I. INTRODUCTION

Achieving fast, automated 3D complex tasks in micromanipulation and micro-assembly with very high accuracy is currently strongly desired in many application fields, such as advanced integrated optics manufacturing. This would indeed make it possible to achieve commercial success with reliable and cost-effective assembly, and packaging technologies. Many emerging industrial on-demand applications made the field of integrated optics particularly challenging, because of their increasing functions and complexity. These applications require the integration of various building components fabricated from different technologies in order to achieve a wide rang of microsystems. The possibility to achieve hybrid nano-photonic components by assembly feeds many original applications. It also brings very relevant solutions to overcome the monolithic integration difficulties as well as to achieve 3D multi functions optical devices manufactured with different technologies [START_REF] Agnus | Robotic microassembly and micromanipulation at femto-st[END_REF], [START_REF] Zimmermann | Packaging and assembly for integrated photonics-a review of the epixpack photonics packaging platform[END_REF], [START_REF] Benson | Assembly of hybrid photonic architectures from nanophotonic constituents[END_REF].

Micro-assembly also aims at optimizing the optical performances. Moreover, it enables the development of complex integrated optical devices, from the assembly of different optical components that then form an integrated optical system. H. Bettahar (e-mails: houari.bettahar@femto-st.fr, olivier.lehmann@femto-st.fr, cclevy@femto-st.fr, philippe.lutz@femto-st.fr, nadege.courjal@femto-st.fr)

Positioning and assembly for integrated optics industry are generally done manually based on passive and active positioning [START_REF] Takigawa | Passive alignment and mounting of linbo 3 waveguide chips on si substrates by low-temperature solid-state bonding of au[END_REF] [5] [START_REF] Sassen | An improved in-plane thermal folded v-beam actuator for optical fibre alignment[END_REF]. In passive positioning, operator may use geometrical devices such as V-grooves [START_REF] Takigawa | Passive alignment and mounting of linbo 3 waveguide chips on si substrates by low-temperature solid-state bonding of au[END_REF]. This kind of approach enables to reach a positioning accuracy of about 1 µm [START_REF] Li | Optimizing fiber coupling with a quasipassive microoptical bench[END_REF], it is accurate enough for many applications as demonstrated for passive alignment of two photonic chips on a silicon optical bench [START_REF] Van Gurp | Passive photonic alignment with submicrometer repeatability and accuracy[END_REF]. There are numerous very accurate machines like flip-chip [START_REF] Yang | A self-aligning flip-chip assembly method using sacrificial positive self-alignment structures[END_REF] [10] that are used in integrated optics industry. However, these technologies are limited to a specific type of products. For complex 3D tasks, positioning based on visual control is usually achieved by teleoperated way in air [START_REF] Kudryavtsev | Stereovision-based control for automated moems assembly[END_REF]. These approaches are more generic, a positioning accuracy of about 10 µm has been achieved [START_REF] Clévy | Positioning accuracy characterization of assembled microscale components for micro-optical benches[END_REF].

However, all these cited strategies rely on geometrical features to relatively position two components relative to one another, the optical performances are not guaranteed and the positioning accuracy is not sufficient. For this sake, the active approach is the mainly used one for integrated optics industry, since it maximizes the optical irradiance during the positioning task. Micro and nano-positioning systems are used to adjust the position of the optical component manually or in a teleoperated way [START_REF] Carminati | Fiber-to-waveguide alignment assisted by a transparent integrated light monitor[END_REF] in order to obtain for instance, the maximum optical irradiance. However, it requires dexterity and skills from the operator along with a lot of precision and attention as it is a tiring task for the operator. The operator applies a non-referenced motion particularly for multi-DOF motion which leads to a local maximum problem on the obtained optical irradiance. The repeatability and accuracy of the positioning may reach some hundreds of nm but at a high cost due to the training of the human operator and the long time needed to achieve the task. Consequently, the performances and the quality of the assembled products are strongly operator-dependent.

Micro-robotic approaches enable to position a frame of a robot, typically an end-effector to the desired location through referenced robotic motions with respect to already attributed reference frames. The extrinsic geometric parameters are the parameters relating the robot end-effector to tool frame, world frame and measurement system frame. The robot end-effector position control is achieved by using the Direct or the Inverse Geometric Model (DGM or IGM) and the transformation matrices between frames. To successfully perform tasks at micro-scale, micro-systems or micro-robots with high positioning accuracy are required. However, this requirement is difficult to obtain at the micro-scale because of many extrinsic geometric parameters acting on micro-robots [START_REF] Shiakolas | On the accuracy, repeatability, and degree of influence of kinematics parameters for industrial robots[END_REF], [START_REF] Veitschegger | Robot accuracy analysis based on kinematics[END_REF]. A solution would be to use additional sensors to provide an outer-loop control (the often used inner-loop feedback one is based on internal sensors of the stages) [START_REF] Clark | Laser-based sensing, measurement, and misalignment control of coupled linear and angular motion for ultrahigh precision movement[END_REF]. However, a key limitation is the complexity of integrating sensors able to measure a relative position with required range, resolution, bandwidth and number of DOF (Degrees-Of-Freedom) [START_REF] Clévy | Signal measurement and estimation techniques for micro and nanotechnology[END_REF].

For this sake, works intend to overcome the active positioning drawbacks by exploiting the high accurate 1-D interferometric measurements that are particularly and naturally produced during classical active positioning, called Fabry-Perot interferences. The aim is then to associate these measurements with a multi-DOF referenced robotic motion. The core objective of this proposed approach is to achieve high accurate multi-DOF positioning in order to get maximum optical performances. In [START_REF] Bettahar | Photorobotic positioning for integrated optics[END_REF], a very preliminary study was conducted considering two optical building components (fiber ferrule, optical component). The study demonstrated that identifying only two parameters (misalignment angles) can largely improve optical performances. This work also stated the potential of a new approach named Photo-robotic because it merges and then correlates the optical, mechanical and robotic behaviors. This work also showed the need to establish a theoretical model and an associated measurement and identification strategy in a formal and generic way to consider more and different types of parameters. For these reasons, this paper aims at generalizing the photo-robotic approach to calibrate all extrinsic geometric parameters of the robot. A new photo-robotic calibration model and an associated measurement strategy are especially proposed. It combines a 1-D Fabry-Perot interferometric measurement with the motions of a multi-DOF robot. Once the robot extrinsic parameters identified, the optical lamella, which is attached to the robot, can be positioned in 6-DOF with very high accuracy. The paper is organized as follows. Section II explains the Fabry-Perot interferometry principle. Section III models geometrically the 6-DOF robot used for positioning. Section IV presents the proposed photo-robotic calibration model. Section V presents the experimental procedure for robot extrinsic parameters calibration. Section VI presents the experimental setup. Section VII investigates the behavior of the calibrated extrinsic parameters of the robot. Section VIII concludes the paper.

II. FABRY-PEROT INTERFEROMETRY

Interferometer sensors are very attracting to use due to their high resolution, accuracy, wide measuring range and bandwidth. These interferometer sensors can be operating based on different interferometry principles such as Michelson and Fabry-Perot interferometry. Active alignment of the optical lamella with respect to the optical fiber produces the Fabry-Perot interferometry phenomenon. The Fabry-Perot interferometry results from the phenomenon of multiple beam interferences that arises when light shines through a cavity bounded by two reflective surfaces as shown in Figure 1. In our case, the involved surfaces are the fiber extremity (R 1 ) in one hand and the lamella surface (R 2 ) in the other hand. Each time the light encounters one of the surfaces, a portion of it is transmitted out, and the remaining part is reflected back. The net effect is to split a single beam into multiple beams which interfere with each other. The system behaves like the Fabry-Perot interferometer, where interferences result from the recombination of multiple reflections between the reflective surfaces R 1 and R 2 respectively.

Based on Fabry-Perot interferometry principle, when the reflective surface (R 2 ) moves along the optical axis (Z-axis), the Fabry-Perot cavity length (L) changes and modifies the reflected light signal accordingly. Therefore, there is a correlation between the FP cavity length L and the interference figure. A slight variation of the FP cavity length can be identified on the interference figure.

The periodicity of the irradiance signal directly depends on the wavelength of the laser signal (λ = 1560 nm). If p is a natural number, the distance between two resonances can be written as in equation 1 (see Figure 2), i.e. the distance variation between two peaks (δL) is a multiple of half wavelength λ of the laser irradiance, such as :

δL = λ 2 • p with p ∈ N (1)
Active positioning and assembly of optical components gives the opportunity to measure a high accurate 1-D relative position along the optical axis with very high resolution. This 1-D Fabry-Perot measure is chosen to be used as a feedback information for positioning of the optical lamella with respect to the fiber. From the other hand, active positioning of optical components generally uses different positioning systems in order to multi-DOF motion, each positioning system realizes motion with respect to it own actuation frame. This leads to a non-referenced multi-DOF motion. The DGM is used to control the different positioning systems as one robot in order to calculate the robot end-effector coordinates with respect to a fixed reference frame as a function of the joints coordinates. The DGM modeling of our 6-DOF robot is the objective of the next section.

III. 6-DOF ROBOT GEOMETRIC MODELING

Geometric model deals with the relationship between the joint and operational coordinates (coordinates with respect to the world reference frame), by considering robot geometric parameters. Geometric model enables to establish of the robot DGM (Direct Geometric Model) and IGM (Inverse Geometric Model), which are successively used to control the robot. It usually consists in determining mathematical description that describes the required coordinate transformations based on different approaches for the geometric representation. One of the most known method has been accomplished by Denavit and Hartenberg [START_REF] Denavit | A kinematic notation for low pair mechanisms based on matrices[END_REF]. The method based on homogeneous transformation matrices. The procedure consists of establishing coordinate systems on each joint axis. Each coordinate system is then related to the next through a specific set of parameters in the homogeneous transformation matrices. Afterwards, many researchers mentioned the model singularity problem raised in D-H representation [START_REF] Ibarra | Determination of linkage parameter and pair variable errors in open chain kinematic linkages using a minimal set of pose measurement data[END_REF]. A new notation has been derived from the popular D-H method by Khalil and Kleinfinger [START_REF] Khalil | A new geometric notation for open and closed-loop robots[END_REF] which can be used to describe the openloop robots and the closed-loop robots with a minimum of parameters and without ambiguities or difficulties. Khalil and Kleinfinger representation is one of the most used conventions for the geometric models. However, in the case of nearly parallel neighboring axes, the common normal is poorly defined and the calibration is ill-conditioned. For this case, Hayati introduced an extra rotational parameter β i about the Y ( i -1 axis [START_REF] Hayati | Improving the absolute positioning accuracy of robot manipulators[END_REF].

The modified version by Khalil and Kleinfinger with additional Hayati parameter β i is used to establish the transformation matrix between two successive links based on two main assumptions: the links are perfectly rigid and the joints are ideal in a sense that there is neither backlash nor elasticity. A serial robot is composed of a sequence of n links and n joints where link 0 is the base of the robot and link n is the terminal link. Joint i connects the link i to the link i -1 and its variable is denoted by θ i . A frame R i is attached to each link i with [START_REF] Khalil | Modeling, identification and control of robots[END_REF]:

• The Z i axis is located along the axis of joint i;

• The X i axis is located along the common normal between the Z i and Z i+1 axes. If Z i and Z i+1 axes are parallel or collinear, the choice of X i is not unique: considerations of symmetry or simplicity then allow a rational choice; • The Y i axis is located using the right-hand rule to form

(X i , Y i , Z i ) coordinate system.
The transformation matrix from frame R i-1 to frame R i is expressed using the following four geometric parameters (Figure 3):

• Joint parallelism: β i is the angle around the axis Y i-1 between the axes Z i-1 and Z i-1 ; h 𝑍 𝑖 𝑌 𝑖 𝑋 𝑖 𝑍 𝑖-1 𝑋′ 𝑖-1 𝑑 𝑖 𝛼 𝑖 𝛽 𝑖 Link 𝑖 Link 𝑖 -1 Joint 𝑖 Joint 𝑖 -1 𝑂 𝑖-1 𝑂 𝑖 𝑑 𝑖+1 𝑟 𝑖 𝜃 𝑖 Joint 𝑖 + 1 𝛽 𝑖 𝛼 𝑖 𝑍′ 𝑖-1 𝑍 𝑖 𝑍 𝑖-1 𝑌 𝑖-1
Figure 3. The additional Hayati parameter β i about Y i-1 is employed for nearly parallel axes.

• Link twist: α i is the angle around the axis X i-1 and between the axes Z i-1 and Z i ; • Link length: d i is the distance along the axis X i-1

between the axes Z i-1 and Z i ; • Joint angle:θ i is the angle around the axis Z i between the axes X i-1 and X i ; • Link offset: r i is the distance along the axis Z i between the axes X i-1 and X i .

The variable of the joint i denoted by q i is θ i if i is rotational and r i if i is prismatic. Hence

q i = θ i (1 -ρ i ) + r i ρ i (2)
where

ρ i = 0, for rotational joint 1, for prismatic joint (3)
Two rotations are required to relate Z i-1 to Z i : a rotation

β i about Y i-1 maps Z i-1 to Z i-1 and X i-1 to X i-1 . Hence, X i-1 axis is located along the common normal between the Z i-1 and Z i axes. A second rotation α i about X i-1 maps Z i-1 to Z i . The angle θ i is now from X i-1 toX i about Z i .
The transformation matrix defining the frame R i in the frame R i-1 is given by

i-1 Ti = RotY i-1 (βi).Rot X i-1 (αi).T rans X i-1 (di) .RotZ i (θi).T ransZ i (ri) (4)
DGM can be obtained by simply concatenating the transformations between the fixed frames of the adjacent links [START_REF] Khalil | Modeling, identification and control of robots[END_REF]:

0 Tn(q) = 0 T1(q1). 1 T2(q2) . . . n-2 Tn-1(qn-1). n-1 Tn(qn) (5)
The inverse problem is to calculate the joint coordinates corresponding to a given operational coordinates of the endeffector. When exists, the form that gives all the possible solutions (there is rarely a single solution) is The Inverse Geometric Model (IGM).

Our objective is to calibrate the extrinsic geometric parameters of the robot. Hence, the intrinsic parameters of the robot are assumed nominals. For this sake, the same origin for all the reference frames is chosen, which corresponds to the real

q 1 q 2 q 3 q 4 q 6 Y w Z w X w q 5 Z 2 Y 5 Z 6 X 4 X 3 X 1 X 5 X 6 Z 1 Z 4 Z 3 Y 4 Y 1 Y 6 Z 5 X 2 Y 3 Y 2 Figure 4
. Robot kinematic model.

Table I NOMINAL PARAMETERS OF DGM OF THE ROBOT BASED ON HAYATI

MODIFIED CONVENTION OF KHALIL AND KLEINFINGER.

Link ρ i β i α i d i θ i r i 1 1 0 0 0 0 q 1 2 1 0 -90 0 -90 q 2 3 1 0 -90 0 -90 q 3 4 0 0 0 0 q 4 0 5 0 0 -90 0 q 5 -90 0 6 0 0 -90 0 q 6 0
rotational center of the robot (intersection of joint 4 and 5) when q 1 = q 2 = q 3 = q 4 = q 5 = q 6 = 0. The world (base) frame R w is assigned to be with the same orientations as the first frame when q 1 = 0. The nominal geometric parameters of the robot are then determined. The assignment of frames is shown in Figure 4, which represents the robot geometric model and the geometrical parameters for each link are given in Table I .

IV. ROBOT EXTRINSIC PARAMETERS CALIBRATION BASED ON 1-D FABRY-PEROT RELATIVE DISTANCE MEASURE

The objective of this section is to achieve the robot extrinsic geometric calibration based on relative distance measure. The distance measure is realized based on very high accurate Fabry-Perot interferometry principle measurement as presented in section II. In order to go through the calibration process, the main reference frames are defined as schematically in Figure 5. R f is assigned to the fiber frame, its origin, labeled O f , is defined as the intersection point between the optical axis and the outer surface of the fiber. R e corresponds to the robot end-effector frame, its origin, labeled O e , corresponds to the real rotational center of the robot (intersection of joint 4 and 5). The ---→ O e X e , ---→ O e Y e and ---→ O e Z e are the expected motion directions of the robot. R t corresponds to the tool frame, its origin, labeled O t , is defined as the intersection between the ---→ O e Z e and the optical lamella plane (Y t O t X t ). R w is assigned to the world frame, its origin, labeled O w , corresponds to the real rotational center of the robot (intersection of joint 4 and 5) when all axes are at zero (initial configuration), with the same directions as the first stage axes R 1 .

DOF robot

Fiber

Optical lamella

𝑅 𝑓 𝑅 𝑡 𝑂 𝑒 𝑂 𝑤 Z 𝑤 X 𝑤 Y 𝑤 X 𝑒 Y 𝑒 Z 𝑒 𝑅 𝑒 𝑅 𝑤 Z 𝑓 X 𝑓 Y 𝑓 O 𝑓 𝑃 𝑟 𝑂 𝑡 𝑍 𝑡 𝑌 𝑡 𝑋 𝑡 𝒘 𝑻 𝒇 𝒆 𝑻 𝒕

End-effector

Figure 5. The photo-robotic positioning scheme with frames assignment.

For simplification, the rotation of the fiber around Y w is assumed to be zero, the rotation of the lamella around Z e axis and the translations of the lamella along X e and Y e axis all are supposed to be zero. The vector ---→ O e Z e is supposed to be collinear with the vector ---→ O w Y w . Thus, the homogeneous transformation matrix from the world frame to the fiber frame is given in equation ( 6), and the transformation matrix from the end-effector frame to the tool frame is given in equation ( 7).

w T f = T ransX w (x wf ).T ransY w (y wf ).T ransZ w (z wf ) .RotX w (β wf ) * RotZ w (θ wf ) =    cos(θ wf )
-sin(θ wf ) 0

x wf cos(β wf ) sin(θ wf ) cos(θ wf ) cos(β wf ) -sin(β wf ) y wf sin(β wf ) sin(θ wf ) cos(θ wf ) sin(β wf ) cos(β wf )

z wf 0 0 0 1    (6)
where x wf , y wf and z wf are the translational coordinates of the origin of the fiber frame O f with respect to the world frame, along X w , Y w and Z w axis respectively. β wf is the angle between Y f and Z w around X w axis and θ wf is the angle between X f and X w around Z w axis.

e Tt = T ransZ e (Zet).RotX e (βet).RotY e (θet) =

   cos(θet) 0 sin(θet) 0 sin(θet) sin(βet) cos(βet) -cos(θet) sin(βet) 0 -cos(βet) sin(θet) sin(βet) cos(θet) cos(βet) Zet 0 0 0 1    (7)
where Z et is the translational coordinate of the origin of the optical lamella O t with respect to the end-effector frame, along Z e axis. β et and θ et are the rotational coordinates of X t and Y t axes of the optical lamella with respect to the end-effector frame, around the X e and Y e respectively. P r is the reflection point projection of the laser spot on the optical lamella plane (Y t O t X t ). w P r ( w P rx , w P ry , w P rz ) is the reflection point projection in the world frame R w . t P r ( t P rx , t P ry , 0) is the reflection point in the tool frame R t .

All the required transformations are defined, and they are used for the proposed extrinsic parameters calibration.

The identification step of the proposed extrinsic geometric parameters calibration consists of three steps:

A. Identification of the optical lamella plane orientation and the optical axis orientation with respect to the world frame (β et , θ et , β wf , θ wf )

The parametric form of the equation of optical axis line is defined by a point on the optical axis line O f and its direction vector -→ Z f . The point O f (x wf , y wf , z wf ) is a point with respect to the world frame R w through which the optical axis line passes. The unit vector -→ Z f (a wf , b wf , c wf ) is the direction vector of the optical axis line with respect to the world frame R w . The 3-dimensional parametric equations of the line carrying the vector ---→ O f Z f (optical axis) can be written as in equation. ( 8):

   x = a wf L + x wf y = b wf L + y wf z = c wf L + z wf (8)
where

a 2 wf + b 2 wf + c 2 wf = 1.
The intersection point between the lamella plane (Y t O t X t ) and the line carries the vector ---→ O e Z e is O t (x wt , y wt , z wt ). The optical lamella plane is defined by the point O t (x wt , y wt , z wt ) in the plane with respect to the world frame and the normal unit vector to the plane -→ Z t (a wt , b wt , c wt ). If M (x, y, z) is any other point which belongs to the lamella plane, then the equation of the lamella plane can be written as follows:

---→ OtM • ---→ OtZt = 0 (9)
This yields:

a wt (x -x wt ) + b wt (y -y wt ) + c wt (z -z wt ) = 0 (10)
By considering d wt = -a wt x wt -b wt y wt -c wt z wt , then:

a wt x + b wt y + c wt z + d wt = 0 (11)
where:

a 2 wt + b 2 wt + c 2 wt = 1.
Differentiating the 3 dimensional optical axis line equation ( 8) and the optical lamella plane equation [START_REF] Kudryavtsev | Stereovision-based control for automated moems assembly[END_REF] at the intersection point between the optical axis line and the lamella plane respectively gives :

   dx = a wf dL dy = b wf dL dz = c wf dL (12)
a wt dx + b wt dy + c wt dz + dd wt = 0

Then by replacing the differentiated optical axis line equation [START_REF] Clévy | Positioning accuracy characterization of assembled microscale components for micro-optical benches[END_REF] in the differentiated optical lamella plane equation [START_REF] Carminati | Fiber-to-waveguide alignment assisted by a transparent integrated light monitor[END_REF] and by replacing dd wt , the equation is yielded:

awta wf dL+bwtb wf dL+cwtc wf dL = awtdx+bwtdy+cwtdz (14)
Equation ( 14) is a non-linear equation relating the relative distance along the optical axis with the 3 relative robot translations. The dx, dy and dz are calculated from the Direct Geometric Model (DGM), dL is calculated using Fabry-Perot measure.

The dedicated measurement strategy for identification of the coordinates of the unit vector Z f (a wf , b wf , c wf ) and the coordinates of the unit vector Z t (a wt , b wt , c wt ) is discussed in the next section.

In order to calculate the (β et , θ et , β wf , θ wf ) angles, the mathematical relationship between the identified (a wt , b wt , c wt , a wf , b wf , c wf ) and the (β et , θ et , β wf , θ wf ) angles need to be established.

The columns of a rotation matrix are orthogonal unit vectors along X, Y and Z axes respectively. Hence, the third column of a rotation matrix represents the unit vector along Z-axis. Accordingly, the transformation matrix w T t of the tool frame in the world frame can be written using the unit vector Z t (a wt , b wt , c wt ) with respect to the world frame as follows:

w Tt = w Te. e Tt =    r r awt r r r bwt r r r cwt r 0 0 0 1

   (15)
The label r is assigned to the matrix elements that are not required in the development of our approach.

Similarly the the transformation matrix e T t of the tool frame in the end-effector frame can be written using the unit vector Z t (a et , b et , c et ) with respect to the end-effector frame as follows:

e Tt =    r r aet 0 r r bet 0 r r cet Zet 0 0 0 1    (16)
The transformation matrix w T e which represents the Direct Geometric Model of the robot (DGM) is supposed nominal. Thus, when q 1 = q 2 = q 3 = q 4 = q 5 = q 6 = 0, the yielded w T e is given:

w Te =    1 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 1    (17)
By substituting equation [START_REF] Li | Optimizing fiber coupling with a quasipassive microoptical bench[END_REF] and equation [START_REF] Clévy | Signal measurement and estimation techniques for micro and nanotechnology[END_REF] in equation [START_REF] Veitschegger | Robot accuracy analysis based on kinematics[END_REF], this yields:

By substituting equation ( 16) and equation [START_REF] Clévy | Signal measurement and estimation techniques for micro and nanotechnology[END_REF] in equation ( 15). This yields:

w T t =     r r a et r r r c et r r r -b et r 0 0 0 1     =     r r a wt r r r b wt r r r c wt r 0 0 0 1     (19)
By matching, we deduce: a et = a wt ; b et = -c wt ; c et = b wt (will be used for the 2nd step).

And by matching, we deduce:

   awt = sin θet bwt = cos θet cos βet cwt = cos θet sin βet (20)
Then, the optical lamella plane orientation angles with respect to the robot end-effector frame can be easily calculated as follow:

βet = atan2(cwt, bwt) θet = atan2(awt cos βet, bwt) (21)
where atan2(b, a) is arc tangent of the two variables a and b.

It is similar to calculating the arc tangent of b/a, except that the signs of both arguments are used to determine the quadrant of the result, which lies in the range [π, π] [START_REF] Slabaugh | Computing euler angles from a rotation matrix[END_REF]. With the same way, using equation ( 22), β wf , θ wf can be deduced by matching.

w T f =    cos(θ wf ) -sin(θ wf ) 0 X wf cos(β wf ) sin(θ wf ) cos(θ wf ) cos(β wf ) -sin(β wf ) Y wf sin(β wf ) sin(θ wf ) cos(θ wf ) sin(β wf ) cos(β wf ) Z wf 0 0 0 1    =    r a wf r r r b wf r r r c wf r r 0 0 0 1    (22)
which gives:

   a wf = -sin θ wf b wf = cos θ wf cos β wf c wf = cos θ wf sin β wf (23)
then the optical axis orientation angles with respect to the world frame can be calculated as follow:

β wf = atan2(c wf , b wf ) θ wf = atan2(-a wf cos β wf , b wf ) (24)
B. Identification of the Cartesian coordinates of the two origins O t (0, 0, Z et ) and O f (x wf , y wf , z wf ) of the optical lamella and the fiber respectively with respect to the world frame When q 1 = q 2 = q 3 = q 3 = q 6 = 0, and by applying two axis-by-axis rotations q 4 and q 5 , equation ( 15) can be written

The new optical lamella plane equation after rotation with respect to the world frame and as shown in Figure 6 is given:

a wt x + b wt y + c wt z + d wt = 0, (26)
          
a wt = cos(q 5 )a et -sin(q 5 )c et b wt = cos(q 4 ) sin(q 5 )a et + sin(q 4 )b et -cos(q 4 ) cos(q 5 )c et c wt = sin(q 4 ) sin(q 5 )a et -cos(q 4 )b et + cos(q 5 ) sin(q 4 )c et d wt = a wt Z et sin(q 5 ) -b wt Z et cos(q 4 ) cos(q 5 ) -c wt Z et cos(q 5 ) sin(q 4 ) (27) By taking the two optical lamella plane equations in the initial pose and after rotation with respect to the world frame as shown in Figure 6 and the optical axis line equation in the world frame as given in equation ( 8), the following equation can be obtained. At the same time this optical lamella distance change on the optical axis after rotation can be measured from the Fabry-Perot interferometry measure.

The parameters (a wt , b wt , c wt , a wf , b wf , c wf ) have already been identified. By applying different rotations q 4 and q 5 , we build a system of linear equations, which will permit to calculate the four parameters (Z et , x wf , y wf , z wf ), the dedicated measurement strategy for identification of O t (0, 0, Z et ) and O f (x wf , y wf , z wf ) is discussed in the next section.

C. Calculation of the reflection point projection P r on the plane (Y t O t X t ) with respect to the tool frame and the world frame ( t P rx , t P ry , w P rx , w P ry , w P rz )

The Cartesian coordinates of O t and O f in the world frame have been identified, the orientations of the optical lamella plane and the optical axis line with respect to the world frame have also been identified. Hence, the point of intersection between the optical axis and the plane of the optical lamella in the world frame can be calculated by calculating the distance l between the optical lamella and the fiber using equation ( 29), and then calculating the Cartesian coordinates of the reflection point w P r ( w P rx , w P ry , w P rz ) in the world frame using:

   w Prx = a wf L + x wf w Pry = b wf L + y wf w Prz = c wf L + z wf (32)
The Direct Geometric Model (DGM) of the robot w T e is supposed nominal, therefore when q 1 = q 2 = q 3 = q 3 = q 4 = q 5 = q 6 = 0, the Direct Geometric Model (DGM) transformation matrix can be deduced as shown in the following equation:

w Te =    1 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 1    (33)
By replacing the already deduced vector w P r and the transformation matrices e T t and w T e , the vector t P r can be deduced as follows:

t P r = e T -1 t * w T -1 e * w P r (34)
In this section, the proposed methodology for extrinsic parameters identification is presented. The next section presents the experimental procedure for robot extrinsic parameters calibration.

V. EXPERIMENTAL INVESTIGATIONS AND BEHAVIOR OF THE CALIBRATED ROBOT

The whole experimental procedure of extrinsic parameters calibration includes measurement and identification, parameters integration and validation steps. Each step is presented in detail.

A. Measurement strategy and identification

There are two adopted measurement strategies that enables to identify two set of different complementary parameters: Notice: in this case, the orientation of the optical lamella and the fiber has to stay constant, which means only translations are allowed. The adopted measurement strategy uses different lamella position combinations and translations. For each lamella position along the axis Z 1 of q 1 = (-150, 0, 150) µm, a continuous translation of 300 µm is applied along X w and a translation of 20 µm along Y w axis, this gives 6 equations. For each lamella position along the axis Z 2 of q 2 = (-150, 0, 150) µm, a continuous translation of 300 µm is applied along X w and Z w , this gives other 6 equations. For each lamella position along the axis Z 3 of q 3 = (-150, 0, 150) µm, a continuous translation of 300 µm is applied along Z w and a translation of 20 µm along Y w axis, this gives again other 6 equations. The overall number of non-linear equations is 18 by applying the 3 translations.

For the translation along X w and Z w axes, the sign of dL is required. For this sake, a wavelength scan is used to calculate the distance L between the fiber and the optical lamella before and after each translation along X w and Z w axes. The difference between the measured distance before and after each translation gives the sign of Fabry-Perot measure of the applied translation.

The free spectral range in wavelength (δλ F SR ) of the Fabry-Perot cavity in general is given by:

δλF SR = λ 2 2nL ( 35
)
where λ is the vacuum wavelength of light. For a linear cavity, such as the Fabry-Perot, L is the length of the cavity, n is the medium index within the cavity. And at the same time, the free spectral range in wavelength (δλ F SR ) can be calculated as follows:

δλF SR = δT λstep tstep (36)
where λ step and t step are the wavelength and time steps respectively. δT can be identified from the reflected irradiance as shown in Figure 7 After calculating the free spectral range in wavelength from equation (36), the cavity distance L can be deduced from equation (35). The resolution of these non-linear equations using nonlinear least-squares solver permits to identify the parameters (a wt , b wt , c wt , a wf , b wf , c wf ), which are the orientations of the optical axis line and the optical lamella plane.

2) Measurement strategy for identification of the Cartesian coordinates of the two origins O t (0, 0, Z et ) and O f (x wf , y wf , z wf ) of the optical lamella and the fiber respectively with respect to the world frame :

Nominal robot IGM ( 𝑒 𝑇 𝑡 𝑡 𝑇 𝑃 𝑟 ) -1

Extrinsic parameters integration

Operational coordinates 𝑋 The adopted measurement strategy uses different combinations of poses. A continuous axis-by-axis rotation is applied around the axis Z 4 of q 4 = -4

Real joints coordinates 𝑞 𝑐

• to 4 • for different angles of q 5 = -4 • , -2 • , 0 • , 2 • , 4 • .
Then, a continuous axis-by-axis rotation is applied around the axis Z 5 of q 5 = -4

• to 4 • for different angles of q 4 = -4 • , -2 • , 0 • , 2 • , 4 • .
The rotations are applied based on velocity control. For each continuous rotation, the lamella distance change along the optical axis is measured, each rotation gives 1 equation. Therefore, the final number of equations is 10.

A system of 10 equations with 4 unknowns is built. The system is overdetermined, since there are more equations than unknowns. The method of ordinary least squares can be used to find an approximate solution to overdetermined systems. For the system AX = b, the least squares formula is obtained from the problem min X ||AX -b||, the solution of which can be written as follows: X = (A T A) -1 A T b [START_REF] Wold | The collinearity problem in linear regression. the partial least squares (pls) approach to generalized inverses[END_REF], where T indicates a matrix transpose.

3) Extrinsic parameters integration: The identification of the previously presented parameters leads to build the transformations t T Pr , w T Pr and e T t . The integration scheme of the identified extrinsic parameters is shown in Figure 8. The inverse kinematic model of the robot is supposed nominal.

Our objective is to rotate the optical lamella around X f , Y f and Z f axes with a fixed reflection point P r , and to translate the optical lamella perpendicular to the optical axis Z f . For this sake, in the validation step, a 100 µm radius circular translation of the lamella along X f and Y f is applied and a 3 rotations are applied around X f , Y f and Z f axes (see Figure 5 for frames assignment).

B. Experimental setup

The experimental set up proposed for 6-DOF automated positioning of optical components is schematically depicted in Figure 9. A 6 DOF serial micro-positioning robot, is proposed to position the optical lamella and control its poses. Further is referred as a robot. An other XYZ manual positioner is fixed in front of the 6 DOF robot, it is used to hold the fiber ferrule and pre-position it. The whole system is referred as photo-robotic system.

The photo-robotic system is mounted on an anti-vibration table. The 6-DOF robot is assembled with minimum geometric errors based on one (XYZ) translation stage followed by a bracket, which is made specifically to join the (XYZ) translation stage (q 1 , q 2 , q 3 ) to a set of two gomiometers (q 4 , q 5 ) and a rotary positioner (q 6 ) (see Figure 9). The two gomiometers are used for rotation about X-axis and Z-axis and the rotary positioner for rotation around Y-axis. A passive XYZ prepositioner is mounted after the rotary positioner for coarse adjustment of the rotation center of the whole assembled robot.

XYZ passive prepositioner

The (XYZ) translation stage is a PIMars P-563.3CD equipped with a E-712 controller. The two goniometers are SGO-60.5 and SGO-77.5, which use the stick-slip principle and is controlled by MCS-3D unit. The SGO-60.5 goniometer's center of rotation is 60.5 mm. The SGO-77.5 goniometer's center of rotation is extended to 77.5 mm. Rotary positioner is SmarAct SR-3610-S which uses the piezo principle and is controlled by MCS-3D unit. All the micro-positioning stages are equipped with internal sensors and are closed-loop controlled in joint layers already.

The optical set up consists in a laser light source (wavelength λ = 1560 nm) connected to a circulator and then to the fiber ferrule. This configuration permits to align and assemble the optical lamella to the fiber. The reflected irradiance from optical lamella through the fiber ferrule is the closed loop information for the robotic micro-manipulator poses control. This set up permits the acquisition of irradiance signal from the reflected light from the optical lamella surface, returning back through the fiber ferrule and then through the beam circulator. Fabry-Perot measurement is used to control the relative position between fiber ferrule/optical lamella actively.

The next section studies and investigates the optical performances after having applied robot extrinsic parameters calibration approach.

C. Behavior study of the calibrated robot extrinsic parameters Figure 10 represents axis-by-axis rotations q 4 , q 5 and q 6 around Z 4 , Z 5 and Z 6 axes respectively before calibration, and square translation on the optical lamella plane before calibration. The fringes (peaks) on the optical irradiance are due to Fabry-Perot cavity variation along the optical axis. This proves the existence of some geometric errors, that can be extrinsic and/or intrinsic. The type of optical irradiance in Figure 10 appears in active positioning of optical components, which leads to local maximum problem. As a result, it is a time consuming especially for multi-DOF alignment process to reach maximum optical power. The objective of robot geometric parameters calibration is to realize 6-DOF positioning without Fabry-Perot cavity variation on the optical axis i.e, realizing rotations around X f , Y f and Z f axes and translation along the optical lamella plane without translating along the optical axis (The reflection point on the lamella plane is kept fixed). For this sake, first extrinsic parameters calibration of the 6-DOF robot is realized. The duration of the measurement and identification phase is about 30 minutes. After achieving measurement and identification phases, the identified extrinsic parameters are obtained and shown in Table II. Then, integration of the obtained extrinsic parameters is performed as shown in Figure 8, where the inverse kinematic model of the robot is supposed nominal. A validation test is then realized which aims at quantifying the performances of the robot after calibration of the extrinsic parameters.

Figure 11 represents the obtained results after having integrated extrinsic parameters for, 3 rotations around X f , Y f and Z f , after identified the reflection point projection on the lamella P r . For rotations around X f , Y f , the corresponding optical irradiance and the 3 translational joints are plotted (since they are the only joints which intervene for these two rotations). For rotations around Z f , the corresponding optical irradiance and the 3 translational joints and the 2 rotational joints (joint 4 and 5) are plotted.

We can notice that after having calibrated robot extrinsic parameters, the joints intervene in order to move the rotation center of robot to the reflection point of the laser spot P r and maintain rotation at that point (reflection point P r ). For rotations around X f and Y f , there exist some fringes (oscillations) on the optical irradiance, which are interpreted by the existence of small Fabry-Perot cavity variation during the two rotations. The estimated Fabry-Perot cavity variation for the two rotations around X f and Y f of rotation range of 6 • is 1.56µm.

For a rotation around Z f , the translational and rotational joints intervene in order to maintain the rotation at the reflection point P r . There exist some fringes on the optical irradiance, which are interpreted by the existence of small Fabry-Perot cavity variation during the rotation. When at least one translational joint is saturated, the frequency of the fringes increase, which represent an increase in Fabry-Perot cavity variation during the rotation. The estimated Fabry-Perot cavity variation for the rotation around Z f is 0.858 µm for a rotation range of 45.5 • .

Figure 12 shows the optical irradiance corresponding to circular translations respectively along the optical lamella plane before and after extrinsic parameters calibration. We can notice that the number of fringes decrease considerably after robot extrinsic parameters calibration. Before calibration, a circular translation of a circumference of 628.32 µm produced 12.87 µm Fabry-Perot cavity variation along the optical axis.

After calibration, the Fabry-Perot cavity variation along the optical axis is about 390 nm. Hence, From the yielded optical irradiance from the circular translation, the optimal position for the optical lamella can be identified to achieve the maximum optical irradiance as highlighted in red color.

VI. CONCLUSION

The paper deals with an integrated optics system achieved through highly accurate multi-DOF robotic active positioning. For this sake, two optical components were chosen, a fiber ferrule and an optical lamella. Fabry-Perot interferometry principle was used to provide a rapid and high accurate measurement.

A photo-robotic approach was proposed for robot extrinsic parameters calibration. This approach relates a 1-D relative distance with a multi-DOF robot motion in order to identify extrinsic geometric parameters based on a photo-robotic calibration model. A dedicated measurement strategy was proposed for high measurement accuracy, where the 1-D relative distance measure is achieved based on the highly accurate 1-D Fabry-Perot interferometry measurement.

A photo-robotic platform dedicated for optical components 6-DOF nano-positioning and assembly was proposed. A 6-DOF serial robot is proposed, which is built based on highly accurate nano-positioning stages commonly used for micromanipulation.

Applying the proposed photo-robotic extrinsic parameters calibration reduced considerably Fabry-Perot cavity variation during rotations around X f , Y f and Z f , after identified the reflection point of the laser spot. We have achieved a Fabry-Perot cavity variation of about 1.56 µm for a rotation range of 6 • around X f and Y f and a Fabry-Perot cavity variation of about 0.858 µm for a rotation range of 45.5 • around Z f . Applying circular translation of 100 µm radius along the optical lamella plane after robot extrinsic parameters calibration reduced also considerably Fabry-Perot cavity variation. We have achieved a Fabry-Perot cavity variation of about 390 nm during a circular translations. From the yielded optical irradiance from the circular translation, the optimal position for the optical lamella can be identified to achieve the maximum optical irradiance.

Figure 1 .Figure 2 .

 12 Photodetector

w

  Tt = w Te. e Tt = βet) sin(θet) sin(βet) cos(θet) cos(βet) Zet -sin(βet) sin(θet) -cos(βet) cos(θet) sin(βet)

5 Figure 6 .

 56 Figure 6. Optical lamella before the new rotation and after the new rotation (around q 5 ) .

  ) cos(q4) 0 sin(q5) sin(q4) -cos(q4) cos(q5) sin(q4wt -Zet sin(q5) r r b wt Zet cos(q4) cos(q5) r r c wt Zet cos(q5) sin(q4

a

  wt x + b wt y + c wt z + d wt = 0 a wt x + b wt y + c wt z + d wt = 0 (28) The intersection of the optical axis line with the first optical lamella plane (before rotation) gives: a wt (a wf L+x wf )+b wt (b wf L+y wf )+c wt (c wf L+z wf )+d wt = 0 (29) The intersection of the optical axis line with the first optical lamella plane (after rotation) gives: a wt (a wf L +x wf )+b wt (b wf L +y wf )+c wt (c wf L +z wf )+d wt = 0 (30) The optical lamella distance change dL along the optical axis after rotation can be deduced as shown follows: dL = L -L = (awtx wf + bwty wf + cwtz wf + dwt) awta wf + bwtb wf + cwtc wf -(a wt x wf + b wt y wf + c wt z wf + d wt ) a wt a wf + b wt b wf + c wt c wf (31)

1 )

 1 Measurement strategy for the identification of the coordinates of the unit vector Z f (a wf , b wf , c wf ) and the coordinates of the unit vector Z t (a wt , b wt , c wt ): Based on Fabry-Perot interferometry principle, when the lamella moves along the optical axis, it produces systematically Fabry-Perot interferometry phenomenon and a change of the Fabry-Perot cavity length (L) happens. Therefore, there is a correlation between the FP cavity length L and the interference figure. A slight variation of the FP cavity length can be identified on the interference figure.

𝛿𝑇

Figure 7 .

 7 Figure 7. Wavelength scan where λstep = 0.01µm and tstep = 0.01s.

Figure 8 .

 8 Figure 8. Extrinsic parameters integration.

Figure 9 . 6 -

 96 Figure 9. 6-DOF Photo-robotic experimental platform.

2 )Figure 10 .

 210 Figure 10. Rotations around Z 4 , Z 5 and Z 6 axes and square translation on the optical lamella plane before calibration.

Figure 11 .

 11 Figure 11. 3 rotations around X f , Y f and Z f axes when the reflection point projection is the reflection point Pr, after robot extrinsic parameters calibration.

Figure 12 .

 12 Figure 12. The optical irradiance corresponding to circular translation on the optical lamella plane before and after extrinsic parameters calibration.

Table II IDENTIFIED

 II EXTRINSIC PARAMETERS. Values 0.744 • -0.595 • 0.732 • -0.520 • 1.52 µm 112.40 µm 900.0 µm 113.3 µm 155.2 µm -2.7 µm -111.34 µm

	Parameters	βet	θet	β wf	θ wf	Zet	x wf	y wf	z wf	w Prx

w Pry w Prz

ACKNOWLEDGMENT

These works have been funded by the Franche-Comté region, nano-roboptic project supported by the EUR EIPHI Graduate School (ANR-17-EURE-0002), ANR CEPAGE ANR-16-CE24-0024 and by the French RENATECH network through its FEMTO-ST technological facility.


Patrick Hénaff
email: patrick.henaff@u-cergy.fr

Vincent Scesa

Fethi Ben Ouezdou
email: fethi.benouezdou@uvsq.fr

Olivier Bruneau
email: bruneau@lisv.uvsq.fr.

Real time implementation of CTRNN and BPTT algorithm to learn on-line biped robot balance: experiments on the standing posture

Keywords: Neural control, Learning algorithms, CTRNN, Legged locomotion, Real-time systems, robotics, Biped robot

This paper describes experimental results regarding the real time implementation of continuous time recurrent neural networks (CTRNN) and the dynamic back-propagation through time (BPTT) algorithm for the on-line learning control laws. Experiments are carried out to control the balance of a biped robot prototype in its standing posture. The neural controller is trained to compensate for external perturbations by controlling the torso's joint motions.

Algorithms are embedded in the real time electronic unit of the robot. On-line learning implementations are presented in detail. The results on learning behavior and control performance demonstrate the strength and the efficiency of the proposed approach.

Introduction

Technological developments have enabled us to build robots with morphologies that are inspired by animals or humans. Therefore, the most recent humanoid robots are technologically complex systems, with an extremely high level of mechanical and electronic integration. They are equipped with complete perceptive systems that enable them to interact with human beings and to move in an environment built for human life.

One of the most important difficulties in controlling humanoid robots is maintaining balance during walking or standing. One solution to this problem is to design controllers based on the zero moment point (ZMP) theory [START_REF] Vukobratovic | Zero-moment-point-Thirty Five years of its life[END_REF]. Another method is to design controllers using bio inspired approaches, i.e., with some capabilities of adaptation and training, leading to the acquisition of reflexes.

Using biologically inspired architectures such as neural networks that are able to learn the "correct" control of the robot's equilibrium is a promising approach. For this purpose, several neural controller-based approaches have been proposed in the past. CMAC (cerebellar model articulation controllers) proposed in 1975 by James [START_REF] Albus | A New Approach to Manipulator Control: the Cerebellar Model Articulation Controller (CMAC)[END_REF], are still studied in the control of legged robots. Recent studies deal with their modeling and generalization properties [START_REF] Horvath | Kernel CMAC With Improved Capability[END_REF] or with their connections to other approaches like fuzzy logic [START_REF] Su | Robust and fast learning for fuzzy cerebellar model articulation controllers[END_REF] or computed torque control [START_REF] Lin | Robust Fault-Tolerant Control for a Biped Robot Using a Recurrent Cerebellar Model Articulation Controller[END_REF]. CMAC have been used to control the balance of biped robots [START_REF] Kun | Adaptive Dynamic Balance of A Biped Robot Using Neural Networks[END_REF] or for robust dynamic walking in simulation [START_REF] Lin | Robust Fault-Tolerant Control for a Biped Robot Using a Recurrent Cerebellar Model Articulation Controller[END_REF]) and for biped robot experiments [START_REF] Sabourin | Robustness of the dynamic walk of a biped robot subjected to disturbing external forces by using CMAC neural networks[END_REF].

Recurrent neural networks (i.e., dynamic neural nets) have been extensively studied in the control of complex systems for many years [START_REF] Marcua | Design of fault detection for a hydraulic looper using dynamic neural networks[END_REF][START_REF] Song | Real-time neural-network midcourse guidance[END_REF]. These Artificial Neural Networks have also been used to design stable walking gaits for biped robots [START_REF] Wu | Robust recurrent Neural Network Control of Biped Robot[END_REF]. Several approaches were based on evolutionary synthesis [START_REF] Fukuda | Stabilization Control of Biped Locomotion Robot based Learning with GAs having Self-adaptive Mutation and Recurrent Neural Networks[END_REF], neural oscillators [START_REF] Taga | Self-organized control of bipedal locomotion by neural oscillators in unpredictable environment[END_REF][START_REF] Geng | A Reflexive Neural Network for Dynamic Biped Walking Control[END_REF], and central pattern generators [START_REF] Righetti | Programmable Central Pattern Generators: an application to biped locomotion control[END_REF][START_REF] Nakanishi | Robot Learning from Demonstration[END_REF]. Recently, researchers have used RNN as predictive compensator [START_REF] Mizunoa | Minimum time ship maneuvering method using neural network and nonlinear model predictive compensator[END_REF] or tracking controllers with selfconstructing properties. Self-constructing algorithms are very interesting approaches because they allow for optimizing on-line the neuronal controller architecture in order to insure the best control performance. In [START_REF] Gao | On-line Adaptive Fuzzy Neural Identification and control of a class of MIMO Nonlinear Systems[END_REF], a self-constructing fuzzy neural controller was proposed for the tracking control of a simulated planar robot manipulator with 2 degrees of freedom. In [START_REF] Hsu Chun-Fei | Intelligent position tracking control for LCM drive using stable on-line selfconstructing recurrent neural network controller with bound architecture[END_REF] a simple growing and pruning algorithm applied to a recurrent neural network has been tested in experimentation to control one degree of freedom of a moving table with a linear ceramic motor system.

Many studies on dynamic neural controllers of robots have focused on continuous-time recurrent neural networks (CTRNN) due to their ability to be universal approximators [START_REF] Beer | Parameter space structure of continuous-time recurrent neural networks[END_REF]. CTRNNs have been used for bio-inspired control because of their abilities to reproduce the full qualitative range of nerve cell phenomenology [START_REF] Beer | Parameter space structure of continuous-time recurrent neural networks[END_REF]. They make it possible to show adaptivity properties based on homeostatic plastic mechanisms [START_REF] Williams | Homeostatic plasticity improves signal propagation in continuous-time recurrent neural networks[END_REF][START_REF] Hoinville | Comparative Study of Two Homeostatic Mechanisms in Evolved Neural Controllers for Legged Locomotion[END_REF]. Moreover, they may be a fine model for generating adaptive behavior because they can learn with the back-propagation through time algorithm (BPTT) [START_REF] Pineda | Generalization of back-propagation to recurrent neural networks[END_REF][START_REF] Werbos | Backpropagation through time: what it does and how to do it[END_REF][START_REF] Pearlmutter | Gradient calculation for dynamic recurrent neural networks: a survey[END_REF]D.E. Rumelhart, Hinto & Williams, 1986;Robinson & Fallside,1987). Unfortunately, these algorithms are complex to implement in real time application, especially for the on-line training of real robots because of the shrinking gradient problem of recurrent neural networks (this gradient instability problem was first studied in [START_REF] Scesa | Time Window Width Influence On Dynamic BPTT(h) Learning Algorithm Performances: Experimental Study[END_REF].

There are still only a few recent in-depth studies on real time implementation, especially with the concept of learning the equilibrium reflex for a biped robot on-line. The scientific objective of the work presented in this paper is to perform an in-depth analysis of the real time performance of continuous time recurrent neural network controllers for learning balance reflexes for a biped robot. In particular, robot abilities to learn online with real time constraints are investigated.

Several contributions in this area have already been published (see [START_REF] Scesa | Dynamic recurrent neural network for biped robot equilibrium control: preliminary results[END_REF]. This paper focuses on describing the real time implementation of the learning algorithms embedded into the robot control unit. Some promising new experimental results are also reported. This paper is organized as follows. The second section deals with the learning models based on the CTRNN and BPTT algorithms. The third section presents the fundamental principles of CTRNN and BPTT.

In the fourth section, the experimental biped platform called ROBIAN is described. The real time implementation of the learning algorithm based on CTRNNs and the back-propagation through time algorithm is detailed Section 5. The sixth section, describes the experimental results in two subparts: a feasibility test for the embedded learning algorithm on controlling the ROBIAN biped torso and a new approach of on-line learning of the equilibrium reflex. Finally, the last section presents the conclusions and potential further developments stemming from this work.

Learning equilibrium with CTRNN

It is well known that the human torso attempts to stabilize the whole body when walking [START_REF] Setiawan | Physical interaction between human and a bipedal humanoid robot: Realization of human-follow walking[END_REF][START_REF] Kubica | Feedforward and deterministic fuzzy control of balance and posture during human gait[END_REF][START_REF] Hyon | Compliant Terrain Adaptation for Biped Humanoids Without Measuring Ground Surface and Contact Forces[END_REF]. It is also known that biped 4 robot balance can be considered to be a global behavior because the upper part (the torso) and the lower part (the legs) of the robot interact with and disturb each other (see Figure 1) during walking. This internal interaction can be taken into account in the synthesis of the biped robot controller. For example, J Morimoto et al. (2006) show that a biped robot can stop and walk using simple sinusoidal desired trajectories with their phases adjusted by a coupled oscillator model. During the walking and halting phases, the robot's balance can be affected by perturbations like external forces applied to its body, including the upper part. One way to control the balance of the robot is to drive the movements of the upper part (i.e., the torso) in order to minimize the perturbations (forces and moments) exerted on the lower part of the robot.

To measure this perturbation, the robot is equipped with a six-force sensor fixed between the torso and the lower part (see Figure 1). Then, if the forces and moments measured by the sensor are close to zero, the equilibrium of the robot is assumed to be controlled [START_REF] Mohammed | Emulation of the Dynamic Effects of Human Torso During a Walking Gait[END_REF]. To achieve this goal, the purpose of this experiment is to learn on-line, using a CTRNN, how to control the robot torso in order to reduce the external perturbations measured by the sensor. This learning is carried out with the biped robot in the standing posture when unknown external forces are applied to its torso.

The control of robot balance must take into account the real features of the mechanism and the external phenomena that are not modeled in the simulation (friction, motor properties, attrition, ground slope and passive prosthesis use…). Information on these phenomena is often only available through their effects on the total energy of the system. Hence, taking into account these phenomena into the equilibrium control is a difficult task. To avoid this obstacle, the controller should be able to adapt its behavior in real time, following a cost function that incorporates information on these phenomena. Moreover, in an optimal control approach, time variations must be included in the adaptation process.

To meet these conditions, dynamic recurrent neural networks with back-propagation through a time learning process was chosen as the most appropriate network. The proposed learning control architecture is shown in Figure 2. The outputs of the network are positions 𝑋, 𝑌, 𝑍 that the torso mechanism then has to reach. Thus, at each time step, the net modifies the actuator speed and accelerations are generated.

Consequently, forces between the upper and lower parts of the robot will be produced to compensate for the external perturbations.

Gradient of cost function

External Perturbations

6

The inputs should show the network the current state of the system. They must be enough representatives to compute a correct control response. They consist of the measured components provided by the 6-forces sensor and the positions of the current motors given by incremental encoders.

To carry out the on-line learning of an optimal control, the parameters of the net must be adapted in real time while the system is running. Back-propagation through time (BPTT), with its ability to integrate the error in the network at each instant is an appropriate solution to solve the parameter adaptation problem. In BPTT, the network is first unrolled in time creating a layer per time step. Then, this algorithm back-propagates the output error on these virtual temporal layers as the classical back-propagation algorithm does on existing ones. The result is the computation of the error gradient in the network and its integration in time. Then, the parameters are modified by the gradient descent algorithm, and the response of the net approaches the desired behavior. Instead of the output error, a cost function that represents the aim of the control can be used for the learning stage. The aim is to bring the forces measured by the sensor closer to zero so that the cost is define by the sum of the squared force components measured by the six axes force sensor in the 𝑋, 𝑌 and 𝑍 directions:

𝐶𝑜𝑠𝑡 𝑡 𝐹 𝑡 𝐹 𝑡 𝐹 𝑡 (1)
This cost function is given for the learning algorithm through a gradient computation for each net output:

𝑒 𝑡 (2)
This equation expresses the influence of each command on the cost function at time t. BPTT will use this gradient to minimize the integral of the cost function (Equation (1)) over the learning time window. Thus, the cost function computed to optimize the control has the following expression:

𝐸 𝐹 𝑡 𝐹 𝑡 𝐹 𝑡 𝑑𝑡 (3)
Instead of computing the exact expression of (3), which would require a strict knowledge of the system, it is possible to use a signal that carries the same information. For the torso experiment, it is defined by the following equations:

𝑒 𝑡 𝐹 𝑡 𝑒 𝑡 𝐹 𝑡 𝑒 𝑡 𝐹 𝑡 (4)

Principles of CTRNN and BPTT

Neurons of a dynamic recurrent neural network are driven by the following equation, derived from [START_REF] Beer | Parameter space structure of continuous-time recurrent neural networks[END_REF][START_REF] Pineda | Generalization of back-propagation to recurrent neural networks[END_REF][START_REF] Pearlmutter | Gradient calculation for dynamic recurrent neural networks: a survey[END_REF]: By the progressive (i.e., explicit) Euler approach, for the ∆𝑡 time step, and the discrete propagation equation depends on the following relations [START_REF] Nguyen | Tau Net: A neural network for modeling temporal variability[END_REF]:

𝑇 𝜕𝑦 𝜕𝑡 𝑦 𝑓 𝑥 (
 Input sum for the 𝑗 neuron :

𝑥 𝑡 𝜔 . 𝑦 𝑡 ∆𝑡 𝑏 𝑡 (6)  Output of the 𝑗 neuron : 8 𝑦 𝑡 ∆𝑡 ∆𝑡 𝑇 . 𝑓 𝑥 𝑡 1 ∆𝑡 𝑇 . 𝑦 𝑡 (7)
In this equation, the ∆ term is a scale parameter. Its value, which cannot be equal to zero, lies in the interval ]0:1] and expresses the response speed of the 𝑗 neuron.

By defining a new variable 𝑆 ∆ , Equation ( 7) can be rewritten as follows:

𝑦 𝑡 ∆𝑡 𝑆 . 𝑓 𝑥 𝑡 1 𝑆 . 𝑦 𝑡 (8)

Truncated Dynamic BPTT learning algorithm

The time constant parameter of each unit reflects the dynamic aspect of the net. However, the learning process must be able to teach it. For this purpose, [START_REF] Pearlmutter | Gradient calculation for dynamic recurrent neural networks: a survey[END_REF] proposed, a dynamical BPTT version for which the time constant, the weights and the bias parameters can be adjusted. The only remaining problem is the memory needed for BPTT algorithms. To compute the exact error gradient, the algorithm has to store all of the network states from the beginning to the current time step. To avoid this excessive memory use, a truncated version, in which the states of the net are only stored during a time window that follows the current instant, can be used [START_REF] Williams | A learning algorithm for continually running fully recurrent neural networks[END_REF][START_REF] Williams | An efficient gradient-based algorithm for on-line training of recurrent network trajectories[END_REF]. Therefore, it only computes an time approximation that approaches the exact gradient as the length of the time window increases. Therefore, the learning algorithm used in the experiments will be a truncated dynamic version of the back-propagation through time algorithm, called TDBPTT.

This learning process begins with an error computation at the output of the net. It represents the remaining gap between the current response and the aim of the learning. This error could simply be a squared difference between the desired and computed response. It could also be the result of a more complex computation, e.g., the equilibrium of a biped robot. The objective of the learning process consists of modifying the network parameters (weights, biases and time scale) in order to minimize a desired criterion. In a control application, this criterion would be the gap between the desired state of the system and the actual state. In an identification process, the criterion would be the error between the neural model and the system being taught.

At first, the algorithm computes an error function that corresponds to the criterion to be minimized. This function, 𝐸, is the error integral between 𝑡 and the current time step, 𝑡, expressed as follows:

𝐸 𝑒 𝜏 . 𝑑𝜏 (9)
where 𝑛 is the number of output neurons, 𝑒 𝜏 is the output error of the 𝑗 neuron stored at time 𝜏 and 𝑡 is the beginning of the integration window. The BPTT algorithm carries out gradient based learning. Thus, the parameters are modified by the negation of the error gradient of each parameter:

∆𝜔 𝜂 ∆𝑏 𝜂 ∆𝑇 𝜂 10
where 𝜂 is the learning step. The parameters of the learning process are the learning step and the time window width 𝑡 𝑡 .

Calculation of the gradient components

The back-propagated cost attached to the 𝑗 neuron is defined by

𝑍 𝜏 𝜕𝐸 𝜕𝑦 𝜏 (11)
and 𝑍 is the back-propagated cost attached to the downstream neuron 𝑘 for each components of the gradient.

10

The calculation of the gradient components for each parameter is carried out on the network "unrolled in time". Each component is estimated for the completely unrolled network, but it results in the values of instantaneous components of the gradients, called the back-propagated gradient components. The total value is the integral, over the unrolled network, of the back-propagated values given by Equation ( 14), for the 𝑗 neuron connected to the 𝑘 neuron by the weight 𝜔 .

By the formula of derivative of the composed functions, the back-propagated components of the gradient can be expressed as follows (see Appendix A for details):

 Gradients on weights:

𝜕𝐸

Calculation of back-propagated costs:

The value of the back-propagated cost 𝑍 expresses the influence of the output of the 𝑗 neuron, at a given time on the current cost. Its computation uses the instantaneous error values 𝑒 𝜏 , calculated outside the network by the comparison of the current state of the network and the objective to be reached.

For each neuron at one instant, 𝜏 of the past, the back-propagation cost comes from three sources. Indeed, the back-propagation "goes up" step-by-step through the neurons from the current activities of the network to The integral, present in the total cost expressed by the Equation ( 9), which made it possible to generate an optimal controller, is consequently distributed over time and over the output neurons. Nevertheless, the calculation of the components of the gradient, reintroduces this integrating aspect, thus guaranteeing the convergence towards an optimal solution over time. This explains the presence of the multiplying term, ∆𝑡, in the expression 𝑍 𝜏 that takes into account the temporal integral of the errors coming from outside.

ROBIAN biped

The learning approach is validated on the ROBIAN biped prototype [START_REF] Konno | Design and development of the biped prototype ROBIAN[END_REF]. ROBIAN consists of two different parts: a locomotion system (lower limbs, i.e, legs) and a special 12 torso mechanism (upper part) that is described in the next section. Each leg has a total of seven dofs [START_REF] Sellaouti | Design and control of a 3DOFs parallel actuated mechanism for biped application[END_REF]: three actuated dofs for the hip, one actuated for the knee, two actuated for the ankle and one passive for the foot, providing a flexible foot system (see Figure 5).

A simplified torso for balance control

The human upper body can be modeled with 13 degrees of freedom, which is not an easy to mechanism control. In order to identify a minimal mechanism of the human upper body that can mimic it when walking, a previous analysis of the six wrench components exerted by the upper part of a virtual manikin on the locomotion apparatus was conducted. It led to the identification of two coupling relations between the upper and lower body. Then, the dynamic equivalence between mechanisms allowed us to identify a trunk mechanism with 4 degrees of freedom that is able to reproduce the dynamic effects of the upper limbs during the walking gait (to understand in detail the advantage of using a simplified torso for balance control, see [START_REF] Mohammed | Emulation of the Dynamic Effects of Human Torso During a Walking Gait[END_REF]). ROBIAN's trunk (see figure 6 and7) is an original mechanism with 4 dofs: one rotational (R) and three translational (P) movements. It is important to note that the three elements, C2, C3 and C4, are directly These dynamic equations show that the movements of the three masses can compensate for unknown external forces applied to the trunk. Thus, in using BPTT to learn this compensation, the cost function and the gradient components for the learning stage will use these equations.

Real time implementation

The real time implementation of the neural controller is a necessary objective to conduct an efficient learning control. The general architecture of the robotic application is based on the exchange of information between different systems: a PC user, an industrial electronic control unit and the ROBIAN robot (Fig. 8.).

The

Multi task organization

The neural networks algorithms, written in the C language, are implemented in the server board. A user interface, written in the Visual Basic graphic language, gives access to all the parameters of the algorithms, to the measurements taken on the robot and to the neural network parameters. The entire application is divided into 5 specific tasks, as shown in Figure 8: the "Main" task, which focuses on most of the information exchanged and manages the objects of the shared memory, the "Measure" task, the "Transfer" task, the "Propagation" task and the "Back-propagation" task.

Scheduling of learning tasks Temporal evaluation

Other experiments have shown that the fastest oscillation period that can be applied to the biped robot in its frontal plan is 500 ms. Thus, the computing time required by the ROBIAN control must be limited to 10 ms in order to maintain sufficient control.

However, the mathematical complexities (number of implied operations) associated with the propagation The algorithmic complexity of the learning tasks is calculated involved in the training, where n is the number of neurons in the net, WT is the Time Window width (number of net states stored). 14 neurons and WT =10 are setup values for the experiments described in Section 5.

Thus, even taking into account the commutation time between the tasks, these two sampling periods remain approximately equal to 2 ms and 50 ms (see next subsection). Figure 9 is a snapshot taken by the numerical oscilloscope of the user interface, which is connected to the electronic control unit. It shows the dependence and synchronization between the propagation and back-propagation tasks. The learning procedure based on the BPTT algorithm uses neuron outputs to calculate backpropagated costs and gradients. However, these same values are modified by the propagation procedure carried out in parallel. It is then not possible to use the same set of variables for the two tasks. Thus, two different tables have been set up in the neuronal algorithms to memorize the states of the network.

Synchronization

For each new iteration of learning, it is necessary to transfer the states from the network generated by the propagation task to the learning task. This procedure, which consists of copying one table into another, is placed upstream of the back-propagation. It is called "temporal photograph" because it memorizes the current instantaneous state of the network. Thus, this table contains the previous states of WT that are present in the temporal window.

The learning procedure gathers four principal stages: temporal photograph (A), backpropagated costs Therefore, their executions are directly related to that of the propagation, as shown in Figure 9. The stage (A) must start only when the propagation (P) is finished. Indeed, it is not necessary for the values that make up a state of the network to be modified when copying this state. To avoid this problem, stage (A) initially checks that the propagation is not running. If the propagation is running, then stage (A) waits until it ends before carrying out the copy. The end of the propagation is announced by a binary semaphore, with a value of 0 (propagation stopped) or 1 (propagation in progress).

In the same way, stage (D) must be carried out when the propagation is stopped. Indeed, this stage modifies the parameters of the network. Thus, it is essential that these parameters are not used when their values are updated. Therefore, this stage (D) waits until the next propagation is finished by observing the semaphore value. If the propagation is not running when launching task (D), it will wait until the propagation launches (propagation in progress) and then finishes (propagation stopped). This will ensure that the propagation will not start during the modification of the parameters. The computation times for this example are given by the table in Figure 10. The duration does not include the waiting periods due to the synchronization process. 6 Experimental Results

Settings for experiments

For this experiment, a 6-5-3 network is used: 6 inputs, 5 hidden neurons and 3 outputs (14 neurons). The number of hidden neurons was selected after some testing. At the beginning of the learning, the weight values of the net are randomly chosen. The time constants are selected from a wide enough range that the net behavior is stable, i.e., that there is no oscillation of the control.

Learning the balance control with Dynamic BPTT

For this approach to the robot's equilibrium control, the experiment is only based on the measured wrench between the upper part and the locomotion system given by the 6-component force sensor. The neural torso controller aims to keep these forces and moments close to zero for the X and Z axes and to weight the value for the Y-axis, even when an external perturbation is applied to the torso. The controller only tries to compensate for external perturbations generated manually on the torso to emulate the coupling effects in a walking gait. The experiment consists of following stages:  ROBIAN stays in the standing posture with all its joints controlled.

 Two successive stages are carried out: learning for 60 seconds and then steering.

 During the learning stage (successive short periods), unknown external forces are manually applied to the torso in X and Z directions (for more comprehension, see Fig. 12, Fig. 13).

 The rotational motion of the torso is locked for this experiment.

 The X and Z moment components are not required for the net computation because their effects can be deduced from the component forces due to the coupling relations (Equations ( 16)). Moreover, as the rotational motion is locked, MY data are also unnecessary.

The learning stage consists of applying forces to each axis, successively. Indeed, each mass is able to compensate only for the forces that appear on its axis. The amplitudes of the applied perturbation and their 22 periods are voluntarily variable in order to show the network different kinds of situations. The values of these external forces are in the maximal range corresponding to the limit of the biped toppling.

On-line learning of disturbance compensation

The aim of the learning stage is to learn to compensate for unknown external disturbances. The time evolutions of responses during the learning stage are depicted in Figures 14,15, and 16. The total duration of the learning stage is 120 seconds. This is the time required by the neural controller to learn the decoupling between the axes when external forces are applied to the ROBIAN torso. The learning algorithm starts at t=4.5s. Until t=11.4s, no external forces are applied, so the measured forces are close to zero (FY close to the ROBIAN torso weight value), and the net does not modify its parameters. Next, during the X-axis perturbation period, the net modifies its X output command to generate a stronger and faster response to the external forces until it becomes insensitive to the disturbance (around t=20s). Then (from t=18 s to 22 s), the Z-axis becomes sensitive to the X-perturbation, and the neural controller explores its control space and learns the decoupling between the two axes. Next (after t=34s), on each axis, the mass has to move in the same direction as the force. Fig. 14. Beginning of on-line learning for the X-axis and Z-axis during time (s). The learning algorithm starts at t=4.5s.The neural controller tries to learn the decoupling between the two axes.

Next, from t=42s (Fig. 14) to t=63s (Fig. 15), it progressively modifies this coupling by decreasing the influence of the X force on its Z control law. From t=62 s, it is clear from Fig. 15 that the network can compensate for each pulse of perturbation due to lateral (Z axis) or frontal (X axis) pushes against the robot.

After this convergence period, the learning continued for about 40 seconds. During that time, the previous changes were carried out (with forces in the opposite direction and more switches between X and Z external perturbation forces). At the end, the network emulated a controller that was able to compensate for a portion of the external forces applied (see Fig. 16). The learning algorithm stops at t=120 s; when a perturbation occurs on an axis, it does not compensate. As far as the Y-axis is concerned, the learning process manages to find motions to compensate for external Y forces in the same way as it does for the two other axes. There was no coupling between this axis and the others. 16. End of on-line learning for the X-axis and Z-axis during time (s). Until t=120s, the neural controller learned the decoupling between the two axis. After, the learning algorithm stops at t=120s, when a perturbation occurs on an axis, it is not compensate. and 𝑍 the back-propagated cost attached to the neuron downstream 𝑘, then each components of the gradient can be written as : 𝑍 𝜏 represents the error coming directly from the output neurons. It does not exist if the neuron is hidden. Hence, for the output neurons only, there is 𝑒 𝜏 0

28

 for a hidden neuron : 𝑍 𝜏 0.

 for an output neuron 𝑍 𝜏 0:

The calculation of the instantaneous errors 𝑒 𝜏 constitutes the coupling between the controlled system and the learning algorithm. This calculation depends on the criterion to minimize. Therefore, it is not possible to write a general formula to obtain these values because they are related to the controlled system and to the cost function:

𝑍

Fig. 1 .

 1 Fig. 1. Control of the biped robot's balance: perturbations caused by unknown external forces can be measured with a force sensor placed between upper and lower part.

Fig. 2 .

 2 Fig. 2. Principle of the on-line learning balance control scheme. Grey lines represent the learning process with cost computation and modification of the net. Each joint position is measured with an incremental electronic sensor. Forces between the upper and lower parts are measured with the six-force sensor. All signals are fed back to the inputs of the control boards of the electronic controller unit (see Section 5).

  the 𝑗 neuron inputs 𝑓 activation function (𝑡𝑎𝑛ℎ or sigmoid) 𝑇 time constant of the 𝑖 neuron 3.1 Propagation algorithm The input data provided to the network are propagated to generate output responses. The propagation depends on the intrinsic parameters of the network and the neurons, which are given as follows for the 𝑗 neuron: 𝜔 weight of the connection from neuron 𝑖 to neuron 𝑗 𝑏 bias of the 𝑗 neuron 𝑇 time constant of the 𝑗 neuron.

Fig. 3

 3 Fig. 3 graphically represents the principles of the propagation through time algorithm. At time 𝑡, each neuron (except input neurons) receives the 𝑡 1 outputs of each other neuron.

Figure 4 Fig. 4 .

 44 Figure 4 represents the calculation of these three elements of the backpropagated cost through the 𝑗 neuron.

Fig. 5 .

 5 Fig. 5.The ROBIAN biped: total height is 1.50 m and weight is 30 kg.

  -10 kg connected to element C1 by prismatic joints (P) and the whole system is connected by a rotational joint (R) to element C0.

Fig. 6 .

 6 Fig.6. The ROBIAN trunk is an RPPP mechanism. On the picture, one can see Y and Z masses that follow their prismatic joint (the X mass is hidden). Total weight is 13 kg.

Fig. 7 .

 7 Fig. 7. Modeling of trunk R3P subjected to external forces, F1 and F2

  communicates with the server board via an Ethernet protocol, which exchanges data with the eight client boards via a special local area network. These client boards exchange with the robot through point-to-point wiring. Programs are run on the server board and are built around the use of objects representing the actuators, the sensors of the robot and their parameters. These objects are written in the shared memory of the electronic unit. The control of each motor is regulated locally via the control boards.

Fig 8 .

 8 Fig 8. Control architecture of the ROBIAN platform

  functions and the training of the neural controller are 𝑂 𝑛 and 𝑂 𝑊 𝑛 , respectively, with n being the number of neurons and WT, the width of the time window used. Consequently, if the two functions are computed sequentially, then the number of operations will evolve out of 𝑂 𝑛 𝑊 𝑛 , and the associated computing time will be able to exceed the limit required. For this reason, the propagation and the backpropagation are separated into two dependent and synchronized tasks. Then, according to the ratio of the computing time (related to the width of the time window) of each task for one training, several propagations in the network are carried out.Table II compares the time performance of propagation and backpropagation. The sampling of the propagation task depends on the computing time necessary for the propagation in the network and, consequently, on the number of parameters defined with the network structure. For example, the propagation time is 2 ms through a completely recurring network consisting of 10 neurons (and 3 ms for 14 neurons and time window of 10). The back-propagation, which is computed along the time window, has an execution time that depends on the width of the time window. For a 10 neuron network with a time window memorizing the 20 preceding states, the back-propagation task will require 50 ms (60 ms for 14 neurons and a time window of 10).

Fig. 9 :

 9 Fig. 9: Visualization of activations of propagation and backpropagation tasks computed in real time by the 10 neurons and there is one backpropagation for 10 propagation step states.

  calculation (B), gradient components computation (C), and parameter modification (D). Stages (B) and (D) are based on the values memorized by the temporal photograph and do not interact with those used by the propagation. On the other hand, stages (A) and (D) copy the network states and modify the parameters.

Figure 10 :

 10 Figure 10: On the left, description of synchronization between the propagation task and the backpropagation procedures. On the right, computation times of the different neural procedures for a 10-neuron network.

Figure 11

 11 Figure11details the execution times of the propagation task, and the other tasks necessary for the learning process. The duration of the induced synchronization waits are also included. In this example, the network consists of 10 neurons and the temporal window has a width of 20 memorized states. The propagation period task is 10 ms, whilst that of the learning task is 100 ms, and the sampling rate of the measurement is 1 ms.Stage (A') corresponds to the 3 ms wait required before the beginning of the temporal photograph (A). The stage (D') represents waiting for 5ms before modifying the network parameters. Note that this stage waits until the propagation starts and finishes before allowing the modification of the parameters (D).

Fig. 12 .

 12 Fig.12. External forces are applied manually during successive short moments following the X direction here (left: pushing; right pulling).

Fig. 13 .

 13 Fig.13. External forces are applied manually during successive short moments following the Z direction here (left: pushing, right: pulling).

Fig. 15 .

 15 Fig. 15. Middle of on-line learning for the X-axis and Z-axis during time (s). The neural controller learns the decoupling between the two axes, and then compensates each perturbation.

Fig.

  Fig. 16. End of on-line learning for the X-axis and Z-axis during time (s). Until t=120s, the neural controller

  electronic unit is built with one server board (SH3 RISC, 96 MHz, 16 Mb of RAM, 128 kb of flash memory, 32 kb of Shared memory) and 8 client boards (SH3 RISC, 96 MHz, 8 Mb for RAM, 128 kb for the flash memory). These 8 boards are dedicated to the control of the 16 DC motors and to the measurement acquisition (2 DC motors per board). The user computer is dedicated to supervising the robot. It

TABLE II

 II

		.
	TIME PERFORMANCES OF PROPAGATION AND BACKPROPAGATION
		Propagation	Backpropagation
	Complexity	O(n²)	O(TW . n²)
	Computational Time(n=10,WT=20)	2 ms	50 ms
	Computational Time(n=14,WT=10)	3 ms	60 ms

  𝑍 𝜏 . 𝑓 𝑥 𝜏 𝛥𝑡 . 𝑦 𝜏 𝛥𝑡 . 𝑑𝜏 (A2) 𝑓 is derivative from the sigmoid activation function 𝑍 𝜏 corresponds to the error coming, at the following instants, from the synapses connecting the neurons:

				𝑆 . 𝑍 𝜏 . 𝑓 ′ 𝑥 𝜏 Δ𝑡 . 𝑑𝜏	(A3)
	 Gradients on scale parameters					 Gradients on weights
	𝜕𝐸 𝜕𝑇	𝜕𝐸 𝜕𝑇	𝜏 . 𝑑𝜏					𝜕𝐸 𝜕𝜔	𝜕𝐸 𝜕𝜔	𝜏 . 𝑑𝜏
		𝜕𝐸 𝜕𝑦 𝜏	𝜕𝑦 𝜏 𝜕𝑇 𝜏	. 𝑑𝜏				𝜕𝐸 𝜕𝑦 𝜏	𝜕𝑦 𝜏 𝜕𝜔 𝜏	. 𝑑𝜏
		𝑍 𝜏		𝜕𝑦 𝜏 𝜕𝑇 𝜏	. 𝑑𝜏				𝑍 𝜏	𝜕𝑦 𝜏 𝜕𝜔 𝜏	. 𝑑𝜏
	𝑆 Δt	. 𝑍 𝜏 . 𝑓 𝑥 𝜏 Δ𝑡	𝑦 𝜏 Δ𝑡 . 𝑑𝜏	(A4)	𝑍 𝜏	𝜕𝑦 𝜏 𝜕𝑥 𝜏	𝜕𝑥 𝜏 𝜕𝜔 𝜏	. 𝑑𝜏
		APPENDIX B Backpropagation in a dynamic neuron.	𝑆 .  Gradients on bias
									𝜕𝐸 𝜕𝑏	𝜕𝐸 𝜕𝑏	𝜏 . 𝑑𝜏
	𝑍 𝜏		𝜕𝐶 𝜕𝑦 𝜏 Δ𝑡	𝜕𝑦 𝜏 Δ𝑡 𝜕𝑥 𝜏	𝜕𝑥 𝜏 𝜕𝑦 𝜏	𝜕𝐸 𝜕𝑦 𝜏	𝜕𝑦 𝜏 𝜕𝑏 𝜏	. 𝑑𝜏
				𝑍 𝜏 Δ𝑡 . 𝑆 . 𝑓 𝑥 𝜏 . 𝜔 𝜏	(B1)	𝑍 𝜏	𝜕𝑦 𝜏 𝜕𝑏 𝜏	. 𝑑𝜏
	 Calculation of 𝒁 𝒋 𝟐 𝝉 :								𝑍 𝜏	𝜕𝑦 𝜏 𝜕𝑥 𝜏	𝜕𝑏 𝜏 𝜕𝑥 𝜏	. 𝑑𝜏

 Calculation of 𝒁 𝒋 𝟏 𝝉 :

  𝑍 𝜏 corresponds to the error coming at the following instants by means of the internal dynamics of the The difference between 𝑍 𝜏 and 𝑍 𝜏 , in the error expression coming from the 𝑗 neuron, comes from the terms considered in the equation (4) because 𝑍 𝜏 takes into account the dynamic part expressed by the term, 1 𝑆 𝑦 𝑡 . Finally, by gathering the equations (B.1), (B.2) and (B.3), the following equation of the backpropagated cost is obtained:

				𝜏	𝜕𝐸 𝜕𝑦 𝜏	𝜕 ∑ 𝜕𝑦 𝜏 𝑒 𝜏	∆𝑡	𝜕𝑒 𝜏 𝜕𝑦 𝜏	∆𝑡	(B2)
	 Calculation of 𝒁 𝒋 𝟑 𝝉 :
	neuron itself :
			𝑍 𝜏	𝜕𝐶 𝜕𝑦 𝜏 ∆𝑡	𝜕𝑦 𝜏 ∆𝑡 𝜕𝑦 𝜏	𝑍 𝜏 ∆𝑡 1 𝑆	(B3)
	𝑍 𝜏	𝑍 𝜏	𝑍 𝜏	𝑍 𝜏
			𝑍 𝜏 Δ𝑡 . 𝑆 . 𝑓 𝑥 𝜏 . 𝜔 𝜏	𝜕𝑒 𝜏 𝜕𝑦 𝜏	. Δ𝑡 𝑍 𝜏 Δ𝑡 . 1 𝑆	(B4)

 Calculation of 𝑍 𝝉 :

External disturbance compensation after learning

Once the learning stage is finished, the net's parameters (biases, weights and time constants) are fixed and stored. The steering consists on using them without adding any new modification.

During the neural control period, the same forces are applied, but the 6-force sensor also measures the influence of the motion of the mass. The difference between the force shapes in the two periods represents the compensation brought about by the mass motion. In a previous study [START_REF] Scesa | Dynamic recurrent neural network for biped robot equilibrium control: preliminary results[END_REF], it was shown that neural control brings a decrease of about 50% in the measured forces on each axis.

Conclusion

This paper addresses the real time implementation of on-line learning control with recurrent Neural Networks. The aim of the control is to compensate for external perturbations due to lateral or frontal pushes against the robot. An experiment is carried out to control the balance of a biped robot in its standing posture.

CTRNN neural nets and BPTT learning have been used to optimize a cost function which explains the rejection of perturbations. The algorithm is embedded in the real time electronic unit of the robot, and the online learning implementation is very detailed. Preliminary results on the learning behavior and the control performances are presented. The implementation enables to focus on the real performance of the neural network controller algorithms and shows that CTRNN are able to learn the balance reflex on-line, when controlling the robot. These experimental results show the strength and efficiency of the proposed approach.

ACKNOWLEDGMENT

The authors would like to thank the French BIA Company, which provides the real time electronic unit that controls ROBIAN and funded the PhD thesis of V. Scesa via a CIFRE contract.

APPENDIX A

Details of gradient components calculation

The backpropagated cost attached to the 𝑗 neuron is defined by: 𝑍 𝜏 𝜕𝐸 𝜕𝑦 𝜏 (A1)