F Aljammal

G Gaborit

Guillaume Revillod

Sylvain Iseni

Lionel Duvillaret

Optical Probe for the Real Time and Vectorial Analysis of the Electric Field Induced by Ionized Gazes

Keywords: optical sensors, plasma diagnostic, intense field, electro-optics, electric field

We here present the potentialities of the electro-optic technique dedicated to the electric field characterization of plasmas. A fully dielectric and millimetre sized probe allows to measure each component of the electric field vector in real time from quasi DC up to several GHz.

1.Introduction

Among all the relevant characteristics of a plasma, its associated electric (E) field behaviour is one of the critical parameter to be analysed. Current and voltage driven by the plasma source are usually measured [START_REF] Schutze | The atmospheric-pressure plasma jet: a review and comparison to other plasma sources[END_REF], but are definitely not sufficient for an exhaustive analysis. Historically, charge density gives additional informations on the plasma and can be assessed thanks to Langmuir probe [START_REF] Hopwood | Langmuir probe measurements of a radio frequency induction plasma[END_REF] or a plasma ion analyser. The non stationary magnetic field can be measured using a B-dot probe [START_REF] Piejak | The electric field and current density in a low-pressure inductive discharge measured with different B-dot probes[END_REF].

Concerning the E-field, numerical simulations [START_REF] Vendor | Numerical modeling of low-pressure RF plasmas[END_REF][START_REF] Mason | An electromagnetic field algorithm for 2D implicit plasma simulation[END_REF] are mainly performed. Actually, there is a lack of available tool allowing the actual measurement in the vicinity or within the plasma.

We here demonstrate the ability of an electrooptic (EO) transducer to provide a complete analysis of the E-field vector induced by the voltage source initiating the ionization process as well as by the plasma itself. The performances of such a technique are in a complete agreement with the requirements for the plasma analysis. Indeed, the measurement bandwidth exceeds several GHz for real time measurement and can be extended up to the THz domain in equivalent time configuration [START_REF] Jianming | Detection of broadband terahertz waves with a laser-induced plasma in gases[END_REF][START_REF] Sanjuan | Full electrooptic terahertz time-domain spectrometer for polarimetric studies[END_REF][START_REF] Gaborit | Emitting and Receiving Terahertz Vectorial Antennas Based on Cubic Electro-Optic Crystals[END_REF]. The dynamic range is greater than 120 dB and ranges from less than 1 V/m to more than 10 6 V/m [START_REF] Gaborit | A nonperturbative electrooptic sensor for in situ electric discharge characterization[END_REF][START_REF] Gaborit | Single shot and vectorial characterization of intense electric field in various environments with pigtailed electrooptic probe[END_REF]. The transverse spatial resolution is better than 1 mm 2 . Finally, the probe includes no metal and is millimetre sized, thus leading to very weak perturbation on the field to be measured, so on the nominal behaviour of the plasma.

2.Electrooptic probe and setup

The optical probe is based on a non centrosymmetric crystal, acting as a linear transducer between the field vector components and its eigen refractive indices. A low noise laser beam (DFB, =1.5 µm) crosses the crystal forth and back and carries, via its polarization state, the information related to the E-field induced modification of the crystal indices. Polarization maintaining fibres allow the remote optical feeding of the EO probe which is connected to an optoelectronic unit. This latter one includes the laser source and also photodiodes ensuring the final conversion into a electrical signal, still proportional to the considered E-field vector component. The probe is packaged with a dielectric sheath to ensure the mechanical and optical integrity. Moreover, this sheath is multi-layered, to improve the adiabatic transition of the field between the outer media and the crystal transducer. The whole measurement scheme is transportable and the performances are not affected by the temperature variation. The probe can be mounted on translation stage to perform the E-field mapping. An example of probes implementation is illustrated in Fig. 1.

3.Experimental results

As an example, we here report some results concerning the analysis of dielectric barrier discharges (DBD). For that purpose, a needle is fed by a 50 Hz high voltage transformer (0.2-25 kVrms). The barriers consist in glass tube and glass plate for the anode and the cathode, respectively. Two probes are used to measure simultaneously the field in the vicinity of both needle and DBD. The obtained temporal waveforms of the radial (horizontal on Fig. 1) field are shown in Fig. 2. The behaviours of the field is analysed for two different feeding voltage. All the measured radial field values are in the 100 kV/m-1.5 MV/m range. When the applied voltage does not induce DBD (Fig. 1a), the measurement in location of probe 1 shows a weaker signal than the probe 2. Indeed, the metallic cylindrical shape of the needle upper part (blue curve) leads to a field aligned with the radial axis; and in between the needle and the ground plate (red curve), the field is mainly vertical. As the measurement is performed and acquired in real time, one can analyse the phase shift  between the two curves.  is here 90° in agreement with a typical RC circuit. When the applied voltage is increased over the threshold (> 16.5 kVrms, Fig. 2b), the radial field remains linear (blue curve) with the applied voltage only for location 2 (close to the metal and far from the discharges). Nevertheless, the red curve presents a typical signature of the discharge once the threshold value is reached, both for positive and negative field. The corresponding peak values exceed 1 MV/m. Moreover the nominal field at 50 Hz has been increased by a factor 5 while the voltage is only doubled. This latter observation is due to modifications of the field direction induced by the DBD, towards a quasi totally horizontal field, along the sensitivity axis of the probe. Concerning the phase shift, the two signals are now in phase because the DBD transforms the previous RC configuration into a pure resistive electrical circuit.

The ability of the probe to measure each component of the field vector independently allows to analyse its temporal evolution in magnitude and orientation. The synchronized acquisition of the radial and longitudinal component, Er and El, leads to polarization trace on the field vector. Fig. 3 gives an example of the vectorial trace of the field in location 1, for two different values of the applied voltage.

4.Conclusion

The performances of the EO probe allow the vectorial and real time analysis of the E-field induced by plasma sources. The properties of the field can be extracted in linear regime as well as for non-linear behaviour. All the properties of the probe will be explained during the conference: sensitivity, linearity, dynamics, vectorial selectivity, frequency bandwidth and so on. Numerous experimental analysis will be presented during the conference: DBD analysis, plasma jet characterization, influence of biological media on the operation of plasma source, ...

Fig. 1 .

 1 Fig. 1. Photography showing two EO probe located in the vicinity of dielectric barrier discharge reactor. Probe 1 analyses the discharges and probe 2 controls the radial field surrounding a high voltage fed needle. The grounded plane is at the bottom of the picture.

  Fig. 2. a) Temporal evolution of the radial E-field in the vicinity of the needle (blue curve) and in between the needle and the ground plane (red curve). The feeding voltage is 10 kV, below the threshold voltage. b) Temporal evolution of the radial E-field in the vicinity of the needle (blue curve) and very close to the DBD. the needle and the ground plane (red curve). The feeding voltage is 20 kV, above the threshold voltage. The inset is a picture of the DBD.

Fig. 3 .

 3 Fig. 3. a) Temporal trace of the field vector deduced from the measurement of Er and El, for an applied voltage of 16 kV, just below the threshold. The ellipse is the feeding curve. b) Temporal trace of the field vector deduced from the measurement of Er and El, for an applied voltage of 20 kV. The ellipse corresponds to the previous one scaled by a factor 20/16=1.25. The expected signature in absence of DBD is an ellipse induced by the contribution of the radial field in zone 2 and 90° phase shifted of the vertical field in the capacitive zone. On can notice a dramatic difference between the signature with and without DBD. In the presence of the DBD, the shape is defined by several almost straight lines. This illustrates a proportional link between Er and El, with linear coefficients depending on the temporal position in the waveform. The left upper zone is clearly associated to the discharge with a huge enhancement of the radial field up to more than 1 MV/m. For lower values of the field strength (center of Fig.3b), Er increases linearly together with El. The slope |Er/El| ranges from 0.62 to 1.16.


Boxiong Zeng

Ting-Li Yang
email: yangtl@126.com

Huiping Shen

Damien Chablat
email: damien.chablat@cnrs.fr

TOPOLOGY DESIGN AND ANALYSIS OF A NOVEL 3-TRANSLATIONAL PARALLEL MECHANISM WITH ANALYTICAL DIRECT POSITION SOLUTIONS AND PARTIAL MOTION DECOUPLING

According to the topological design theory and method of parallel mechanism (PM) based on position and orientation characteristic (POC) equations, this paper design a novel 3translation (3T) PM that has three advantages, i.e., ① it consists on three actuated prismatic joints, ② the PM has analytical direct position solutions, and ③ the PM is of partial motion decoupling property. Firstly, the main topological characteristics such as the POC, degree of freedom and coupling degree are calculated for kinematics modelling. Due to the special constraint feature of the 3-translation, the analytical direct position solutions of the PM can be directly obtained without needing to use one-dimensional search method. Further, the conditions of the singular configuration of the PM, as well as the singularity location inside the workspace are analyzed according to the inverse kinematics.

INTRODUCTION

In many industrial production lines, process operations require pure translation movements only. Therefore, the 3-DOF translational parallel mechanism (TPM) has a significant potential value due to its small number of actuated components, relatively simple structure and easily to be controlled.

Many scholars have being studied the TPM. For example, original design of 3-DOF TPM is the Delta Robot which was presented by Clavel [START_REF] Clavel | A Fast Robot with Parallel Geometry[END_REF]. The Delta-based structure manipulators of TPM have been developed [START_REF] Stock | Optimal Kinematic Design of Spatial Parallel Manipulators: Application to Linear Delta Robot[END_REF][START_REF] Bouri | The Linear Delta: Developments and Applications[END_REF][START_REF] Kelaiaia | Multiobjective optimization of a linear Delta parallel robot[END_REF]. Tsai et al [START_REF] Tsai | Kinematics of a Novel Three DoF Translational Platform[END_REF][START_REF] Tsai | Kinematics and optimization of a spatial 3-UPU parallel manipulator[END_REF] presented the 3-DOF TPM, the moving actuators of which are prismatic joints and the sub-chain is 4R parallelogram mechanism (P is prismatic joint and R is revolute joint). In [START_REF] Zhao | Kinematics analysis of a three dimensional mobile parallel platform mechanism[END_REF][START_REF] Yin | Workspace Analysis of 3-DOF translational 3-RRC Parallel Mechanism[END_REF], the authors suggested a 3-RRC TPM and developed the kinematics and workspace analysis (C is cylindrical joint). Kong et al [START_REF] Kong | Kinematics and singularity analysis of a novel type of 3-CRR 3-DOFtranslational parallel manipulator[END_REF] proposed a 3-CRR mechanism with good motion performance and no obvious singular position. Li et al [START_REF] Li | Kinematics of a Special 3-DOF 3-UPU Parallel Manipulator[END_REF][START_REF] Li | Kinematic characteristics of a special 3-UPU parallel platform manipulator[END_REF] developed a 3-UPU PM (U is universal joint) and analyzed the instantaneous kinematics performance of the TPM. Yu et al [START_REF] Yu | Comprehensive Research on 3-DOF Translational Parallel Mechanism[END_REF] carried out a comprehensive analysis of the three-dimensional TPM configuration based on the screw theory. Lu et al [START_REF] Lu | Motion and Workspace Analysis of a Novel 3-Translational Parallel Mechanism[END_REF] proposed a 3-RRRP (4R) three-translation PM and analyzed the kinematics and workspace. Yang et al [START_REF] Yang | Topology Structure Design of Robot Mechanisms[END_REF][START_REF] Yang | Theory and Application of Robot Mechanism Topology[END_REF][START_REF] Yang | Topology Design of Robot Mechanism[END_REF][START_REF] Yang | Composition Principle Based on Single-Open-Chain Unit for General Spatial Mechanisms and Its Application[END_REF] based on the single opened chains (SOC) units theory to synthesize the 3T0R PM, a variety of new TPMs were synthesized and then classified, Considering the anisotropy of kinematics, Zhao et al [START_REF] Zhao | Dimensional synthesis of a three translational degrees of freedom parallel robot while considering kinematic anisotropic property[END_REF] analyzed the dimension synthesis and kinematics of the 3-DOF TPM based on the Delta PM. Zeng et al [START_REF] Zeng | Tri-pyramid Robot: Design and kinematic analysis of a 3-DOF translational parallel manipulator[END_REF][START_REF] Zeng | Tri-pyramid Robot: stiffness modeling of a 3-DOF translational parallel manipulator[END_REF][START_REF] Lee | Error modeling for sensitivity analysis and calibration of the tri-pyramid parallel robot[END_REF] introduced a 3-DOF TPM called as Tri-pyramid robot and presented a more detailed analytical approach for the Jacobi matrix. Prause et al [START_REF] Prause | Comparison of Parallel Kinematic Machines with Three Translational Degrees of Freedom and Linear Actuation[END_REF] compared the characteristics of dimensional synthesis, boundary conditions and workspace for various 3-DOF TPM for the better performance. Mahmood et al [START_REF] Mahmood | Kinematic Analysis and Design of a 3-DOF Translational Parallel Robot[END_REF] proposed a 3-DOF 3-[P2(US) mechanism and analyzed its kinematics and dexterity.

However, the most previous TPMs generally suffer from two major problems: i) the degree of coupling κ of these PMs is greater than zero, which means its direct position solution is generally not analytical, and ii) these PMs do not have inputoutput decoupling characteristics [START_REF] Shen | Kinematic decoupling design method and application of parallel mechanism[END_REF], which leads to the complexity of motion control and path planning.

According to topology design theory of PM based on position and orientation characteristics (POC) equations [START_REF] Yang | Theory and Application of Robot Mechanism Topology[END_REF][START_REF] Zhao | Dimensional synthesis of a three translational degrees of freedom parallel robot while considering kinematic anisotropic property[END_REF], a new TPM is proposed in this paper. The TPM is designed with simple structure and features one coupling degree, it consists of prismatic and revolute joints, and has analytical direct position solutions and partial motion decoupling property. The position solutions, singularity, workspace and its internal singularity of the PM are analyzed.

DESIGN AND TOPOLOGY ANALYSIS

Topological design

The 3T parallel manipulator proposed in this paper is illustrated in Fig. 1. The base platform 0 is connected to the moving platform 1 by two hybrid chains that contain looped Simple-Opened-Chains. Such hybrid chains are called HSOCs, their structural and geometric constraints are as follows:

FIGURE 1 KINEMATIC STRUCTURE OF THE 3T PM

(1) For the 6-bar planar mechanism loop (abbreviation: 2P4R planar mechanism) in right side of Fig. 1, two revolute joints R 3 and R 4 whose axes are parallel to each other are connected in series, where R 3 is connected to link 11 and R 4 is connected to the moving platform 1 to obtain the first HSOC branch (denoted as: hybrid chain I). Two prismatic joints P 1 and P 2 will be used to be actuated.

(2) The left side branch is made up of a prismatic joint P 3 and two 4R parallelogram mechanisms connected in series, and the parallelograms connected from P 3 to the moving platform 1 are respectively recorded as ①, ②, after P 3 and the parallelogram ① are rigidly connected in the same plane, they are connected to the parallelogram ② in their orthogonal plane to obtain the second HSOC branch (denoted as: hybrid chainⅡ).

(3) The prismatic joints P 1 , P 2 and P 3 are connected to the base platform 0; P 1 and P 2 are arranged coaxially, and prismatic P 1 is parallel to P 2 . When the PM is moving, the 2P4R planar mechanism is always parallel to the plane of the parallelogram ①.

Analysis of topology characteristics

Analysis of the POC set

The POC set equations for serial and parallel mechanisms are expressed respectively as follows:

 m i Ji bi M M 1   ( 1 ) bi n i Pa M M  1   ( 2 )
where Ji M -POC set generated by the i-th joint. bi M -POC set generated by the end link of i-th branched chain. Pa M -POC set generated by the moving platform of PM.

Obviously, the output motion of the intermediate link 11 of the 2P4R planar mechanism on the hybrid chain I is two translations and one rotation (2T1R). The output motion of the link S of the parallelogram ① on the hybrid chain II is two translations (2T), and the output motion of the link T of the parallelogram ② is three translations (3T). Therefore, the topological architecture of the hybrid chain I and II of the PM can be equivalently denoted as, respectively:

       4 3 ) 4 2 ( ) 4 2 ( 2 ) 4 2 ( 1 1 || ) , ( R R R P P HSOC R P R P R P       ) 4 ( ) 4 ( 3 2 R R P P P HSOC
The POC sets of the end of the two HSOC S are determined according to Eqs. ( 1) and (2) as follows:

                       )) , ( (|| ) (|| ) ( ) (|| ) ( 3 12 2 3 3 1 3 2 12 1 12 2 1 R R r t R r R t R r R t M HSOC                        0 3 0 1 2 1 0 3 1 )) ( (|| ) (|| 2 r t r R R R R t r P t M di ci bi ai i HSOC 
The POC set of the moving platform of this PM is determined from Eq. (2) by

  2 1 HSOC HSOC Pa M M M        0 3 r t
This formula indicates that the moving platform 1 of the PM produces three translations motion. It is further known that the hybrid chain II in the mechanism itself can realize the design requirement of three translations, which simultaneously constrains the two rotational outputs of the hybrid chain I.

Determining the DOF

The general and full-cycle DOF formula for PMs proposed in author's work [START_REF] Yang | Topology Design of Robot Mechanism[END_REF] is given below:

      m i v j Lj i f F 1 1  ( 3 )            v j b b j i Lj M M j i 1 1 ) 1 ( . dim   ) (  ( 4 )
where

F -DOF of PM. f i -DOF of the i th joint.
m -number of all joints of the PM.

v -number of independent loops, and v=m-n+1.

n -number of links. j L  -number of independent equations of the j th loop.

 j i b i M 1 
-POC set generated by the sub-PM formed by the former j branches.

( 1) bj M  -POC set generated by the end link of j+1 sub-chains.

The PM can be decomposed into two independent loops, and their constraint equations are calculated as follows:

① The first independent loop is consisted of the 2P4R planar mechanism in the hybrid chain I, the LOOP 1 is deduced as:

      ) 4 2 ( ) 4 2 ( 2 ) 4 2 ( 1 1
) , (

       3 ) 4 ( ) 4 ( 4 3 2 || P P P R R LOOP R R
In accordance with Eq.( 4), the independent displacement equation number 2 L  of the second loop can be obtained as follows:

5 )) , ( (|| . dim ) (|| ) (|| ) ( . dim 3 12 2 3 3 1 3 12 1 12 2 2                                            R R r t R r t R r R t L   Thus, the DOF of the PM is calculated from Eq. (3) as            2 1 1 3 ) 5 3 ( ) 5 6 ( j L m i i j f F 
Therefore, the DOF of the PM is 3, and when the prismatic joints P 1 , P 2 and P 3 on the base platform 0 are the actuated joints, the moving platform 1 can realize 3-translational motion outputs.

Determining the coupling degree

According to the composition principle of mechanism based on single-opened-chains (SOC) units, any PM can be decomposed into a series of Assur kinematics chains (AKC), and an AKC with v independent loops can be decomposed into v SOC. The constraint of the j th SOC is defined [START_REF] Yang | Topology Design of Robot Mechanism[END_REF][START_REF] Yang | Composition Principle Based on Single-Open-Chain Unit for General Spatial Mechanisms and Its Application[END_REF] by

                                , 3 , 2 , 1 0 1 , 2 , 4 , 5 0 1 j j j L j m i i j j j I f  ( 5 )
where j m -number of joints contained in the j th SOC j .

i f -DOF of the i th joints. j I -number of actuated joints in the j th SOC j . j L  -number of independent equations of the j th loop.

For an AKC, it must be satisfied

0 1     v j j
then, the coupling degree of AKC [START_REF] Yang | Topology Design of Robot Mechanism[END_REF][START_REF] Yang | Composition Principle Based on Single-Open-Chain Unit for General Spatial Mechanisms and Its Application[END_REF] is

          v j j in 1 m 2 1  ( 6 )
The physical meaning of the coupling degree κ can be explained in this way. The coupling degree κ reflects the correlation and dependence between kinematic variables of each independent loop of the mechanism. It has been proved that the higher κ, the greater the complexity of the kinematic and dynamic solutions of the mechanism will be.

The independent displacement equations of LOOP 1 and LOOP 2 have been calculated in the previous section Determining the DOF, i.e.,

3 1  L  , 5 2  L 
, thus, the degree of constraint of the two independent loops are calculated by Eq. ( 5), respectively, can be obtained as follows:

1 3 2 6 1 1 1 1 1           L m i i I f  1 5 1 5 2 2 2 1 2            L m i i I f



The coupling degrees of the AKC is calculated by Eq. ( 6) as

1 ) 1 1 ( 2 1 2 1 1          v j j k
Thus, the PM contains only one AKC, and its coupling degrees equals to 1. Therefore, when solving the direct position solutions of the PM, it is necessary to set only one virtual variable in the loop whose degree of constraint is positive ) 0 (   j . Then, a position constraint equation with this virtual variable is established in the loop with the negative constraint ) 0 (   j , and the real value of the virtual variable can be obtained by the one-dimensional search method, accordingly to obtain the direct position solutions of the PM.

However, owing to the special 3-translational feature constraint of the PM, the loop with the negative degree of constraint ) 0 (   j can be directly applied to the geometric constraint of the loop with the positive degree of constraint ) 0 (   j , i.e., The motion of the link 11 is always parallel to the base platform 0, from which the virtual variable is easily obtained, and there is no need to solve the virtual variable by one-dimensional search method. Therefore, the analytical direct position solutions of the PM can be directly obtained in the following section, which greatly simplifies the process of the direct solutions. This method of solving the direct solutions of the virtual variable directly by special geometric constraints is of general.

POSITION ANALYSIS

The coordinate system and parameterization

The kinematics modeling of the PM is shown in Fig. 2. The base platform 0 is a rectangle having a length and a width of 2a and 2b respectively. The frame coordinate system O-XYZ is established on the geometric center of the base platform 0, the X and Y axes of which are perpendicular and parallel to the line A 1 A 2 , and the Z axis is determined by the right hand Cartesian coordinate rule. The moving coordinate system O'-X'Y'Z' is established at the center of the moving platform 1, the X' and Y' axes of which are coincide and perpendicular to the

LOOP (PARTIAL) IN THE XOZ DIRECTION FIGURE 2 KINEMATIC MODELING OF THE PM

The angle between the vectors B 1 C 1 and the Y axis is  , and the  assigned as virtual variable. The angles between the vectors D 1 D 2 , D 3 E 3 and the X axis are  and  respectively.

Direct kinematics

To solve the direct kinematics, i.e., to compute the position O'(x,y,z) of the moving platform when setting the position values of the prismatic joints P 1 , P 2 and P 3 (with the

1)Solving the first loop with positive degree of constraint

1 LOOP : 2 2 2 1 1 1 A B C C B A     
The coordinates of points A 1 , A 2 and A 3 on the base platform 0 are respectively

T A y b A ) 0 , , ( 1 1   , T A y b A ) 0 , , (
C C z z  ( 7 )
Therefore, the coordinates of points C 1 and C 2 are calculated as

T A l l l y b C ) sin , cos , ( 2 1 2 1 1       T A l l l l y b C ) sin , cos , ( 2 1 3 2 2 1       
Due to the link length constraints defined by Thus, the second loop with negative degree of constraint acts on the special geometric constraint Eq. ( 7) on the first loop with positive degree of constraint, which is the key to directly finding the analytical solutions of the virtual variable  . This is an advantage for easy obtaining the analytical direct solutions from the topological constraint analysis of the PM.

2 ) Solving the second loop with negative degree of constraint 2 LOOP :

3 3 3 3 3 3 2 1 A B C D E F D D       
The coordinates of points D 1 and D 2 obtained from points C 1 and C 2 are calculated as

T A l l l l y b D ) sin , 2 / cos , ( 2 1 3 2 1 1                             sin sin 2 / cos
) ( 2 2 1    H t H 2 1 H H t    (14)
where , sin

7 8 2 1 l l l H     . ) ( 2 2 6 2 3 3 B C y y l H   
When the PM is moving, the 2P4R planar mechanism is always parallel with the plane of the parallelogram ① , therefore, there is always relation as After eliminate  from Eqs. ( 13) and ( 16), there is

0 cos sin 3 2 1    J J J   let 2 tan 1   k 3 2 2 3 2 2 2 1 1 arctan 2 J J J J J J       (17) where . ) ( 4 ), ( 4 , 2 2 2 2 4 2 6 3 4 2 4 1 d b t l l J d b l J t l J        
Finally, putting the values of  and  obtained by Eqs. ( 9) and ( 17) respectively, into Eq. ( 10), the coordinates of point O'on the moving platform 1 can be obtained.

From Eq. ( 9), ) , ( 2 For the convenience of understanding, the above calculation process can be described by Fig. 3.

                                                                                    ' ' ' 3 3 3 3 2 2

FIGURE 3 FLOW CHART OF DIRECT POSITION SOLUTIONS

It can be seen that the geometric constraint Eqs. ( 7), ( 15) and ( 16) are the key to find the analytical formula of the first and second loop position equations of the PM.

Inverse kinematics

To solve the inverse kinematics, we compute the values of For a given position of the moving platform, from Eqs. [START_REF] Li | Kinematics of a Special 3-DOF 3-UPU Parallel Manipulator[END_REF] and ( 16), the angles  and  are calculated as Further, the coordinates of points C 1 and C 2 are defined as:

T l z l y b C ) sin , 2 / , ( 4 3 1      T l z l y b C ) sin , 2 / , ( 4 3 2     
In addition, the coordinates of point C 3 have been given by Eq. [START_REF] Li | Kinematic characteristics of a special 3-UPU parallel platform manipulator[END_REF]. Therefore, due to the link length constraints defined by [START_REF] Zeng | Tri-pyramid Robot: stiffness modeling of a 3-DOF translational parallel manipulator[END_REF] From Eqs. [START_REF] Zeng | Tri-pyramid Robot: stiffness modeling of a 3-DOF translational parallel manipulator[END_REF], we can evaluate )

                         2 6 2 2 2 2 2 2 2 2 2 2 2 2 2 ) ( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( 3 3 3 3 3 3 2 2 2 2 2 2 1 1 1 1 1 1 l z z y y x x l z z y y x x l z z y y x x B C B C B C B C B C B C B C B C B C
i A as following ) 3 , 2 , 1 (    i M y y i C A i i ( 21) where , ) ( , ) ( 2 1 2 2 2 2 1 2 2 1 2 1 l z l M l z l M C C       . ) ( 3 , 2 , 1 (  i y

Numerical valuation for direct and inverse solutions

Direct solutions

The dimension parameters of the PM are (unit: mm): It can be seen that the inverse solutions data of the No. 4 group in Table 2 is consistent with the three input values given when the direct solution is solved, which proves the correctness of the direct and inverse solutions.

300  a , 150  b , 50  d , 30 1  l , 280 2  l , 1 162

SINGULARITY ANALYSIS

Method of singularity analysis

This paper uses the Jacobian matrix method to analyze the singularity configuration of the PM. Taking the first derivative of time t from Eqs. ( 18) and ( 19), we have

  sin 4 l x     ( 22)   sin 6 l x     ( 23)
Then, taking the first derivative of time t from length constraint Eq. ( 20), and then substituting Eqs. ( 22) and ( 23) into the equations, there are

) 3 , 2 , 1 ( 0 3 2 1          i y u z f y f x f i A ii i i i ( 24)
where . , , Whether the matrices p J and q J are singular, the singularity of the PM is divided into the following three categories:

3 3 2 2 1 1 33 22 11 B C B C B C y y u y y u y y u       ). ( ), ( ), ( cot 1 1 1 1 1 1 13 12 11 B C B C B C z z f y y f z z f        ). ( ), ( ), ( cot 2 2 2 2 2 2 23 22 21 B C B C B C z z f y y f z z f        ). ( ), ( ), ( cot
(1) When 0 ) det(  q J

, the PM has serial singularity .

(

) When 0 ) det(  p J 2
, the PM has parallel singularity.

(

) When 0 ) det( ) det(   p q J J 3
, the PM has comprehensive singularity.

The result of singularity analysis

(1) Serial singularity: When the PM has serial singularity, it means that each of the two branches near the drive link is folded or fully deployed. At this movement, 0 ) det(  q J , the DOF of the moving platform is reduced, and the set W of the equation solution is:

  3 2 1 W W W W    (26) where   0 1 1 1    B C y y W
, that is, the three points A 1 , B 1 and C 1 are collinear.

  0 2 2 2    B C y y W
, that is, the three points A 2 , B 2 and C 2 are collinear.

 

0 3 3 3    B C y y W
, that is, the three points A 3 , B 3 and C 3 are collinear.

Example: The 3D configuration of satisfying condition W 3 is shown in Fig. 5.

FIGURE 5 EXAMPLE OF SERIAL SINGULARITY CONFIGURATION

(2) Parallel singularity: When the PM has parallel singularity, it means that each branch is close to the link of the moving platform in a state of being folded together or fully deployed. At this movement, the DOF of the moving platform is increased, and even if the input link is locked, there may be has DOF output on the moving platform, assuming:

  ) 3 , 2 , 1 ( 3 2 1    i e f f f i i i i
f f f k f f f  is satisfied, there are: ), ( cot ) ( cot 3 3 1 1 B C B C z z k z z      ). ( ) ( 3 3 1 1 B C B C z z k z z   
f f f k f f f k f f f  
The calculation by MATLAB shows that the solution of k cannot be solved in this case, Therefore, this situation does not exist.

(3) Comprehensive singularity:

0 ) det( ) det(   p q J J
, the serial singularity and the parallel singularity occur simultaneously. In this configuration, the actuated joints and the end effector of the PM have non-zero inputs and outputs that do not affect each other instantaneously, and the corresponding pose is the third kind of singularity. In this kind of singularity, the PM will lose the freedom and it should be avoided during the design phase of the PM.

WORKSPACE ANALYSIS

This paper uses the limit boundary search method to analyze the workspace of the PM, that is, the search range of the workspace is first assigned according to the length of the link. Then, based on the inverse position solution, all the points satisfying the constraint are searched, and the 3D map composed of these points is the workspace of the PM. The serial singularity can be avoided by actual control, so this paper mainly discusses the parallel singularity trajectory. According to the link length constraint of Eq. ( 20) and the discriminant 0 ) det(  p J , the parallel singularity trajectory can be obtained, as shown in Fig. 7. Among them, the green part is the non-singularity workspace, and the red part is the singularity area, which indicates that there is a large nonsingularity area inside the workspace. Fig. 8 is the projection view of the workspace in the XOZ and YOZ directions.

CONCLUSIONS

The new 3-translational (3T) PM proposed in this paper has three advantages: (1) it is only composed of three actuated prismatic joints and other passive revolute joints, which is easy to be manufactured and assembled; (2) it has analytical direct position solutions, which brings the great convenience to error analysis, dimensional synthesis, stiffness and dynamics research; and (3) it has partial input-output motion decoupling, which is very beneficial to the trajectory planning and motion control of the PM.

According to the kinematics modeling principle proposed by the author based on the single-opened-chains method, in the first loop with positive constraint, the set one virtual variable  can be directly obtained by the special geometric constraint condition that the output link of the first loop always maintains the horizontal position (the condition is provided by the second loop with negative constraint). Therefore, the entire analytical position solutions are obtained without solving the virtual variable  by the geometric constraint equation in the second loop with negative constraint. This is the advantage of the topology of the PM being different from other PM, and it has analytical direct solutions. The method has clear physical meaning and simple calculation.

Based on the inverse solution, the conditions and locations of the three types of singularity configurations of the PM are obtained, and the size of the workspace of the PM and its parallel singularity area are given. The work of this paper lays the foundation for the stiffness, trajectory planning, motion control, dynamics analysis and prototype design of the PM.

②

  The above 2P4R planar mechanism and the following sub-

line D 2 F

 2 3 , and the Z' axis is determined the right hand Cartesian coordinate rule. The length of the three driving links 2 is l 1 , the length of the connecting links 9 and 10 on the hybrid chain I is l 2 , and the lengths of the intermediate links 11 and 12 are l 3 , l 4 respectively. The length of the parallelogram short links 3, 6 on the hybrid chain II is l 5 . the point B 3 , C 3 , D 3 and E 3 are the midpoint of the short edge, the length of the long links 4, 7 is l 6 , The length of the connecting link 5 between the parallelograms is l 7 . The length of the connecting link 8 is l 8 , and the length of the line D 2 F 3 on the moving platform 1 is 2d.

.

  each end-point of the three same actuated links 2, i.e., B 1 , B 2 and B 3 are easily calculated as Due to the special constraint of the three translations of the moving platform 1, during the movement of the PM, the intermediate link 11 of the 2P4R planar mechanism is always parallel to the base platform 0, that is, C 1 C 2 || A 1 A 2 , then we have 2 1

  can be deduced as below.

  of  cannot be determined at this time, and the PM has parallel singularity, this situation of which should be avoided.

  coordinates of points F 3 , E 3 , D 3 and C 3 are represented by the coordinates of point O' defined as: can be deduced as below.

  ., the PM has partial input-output motion decoupling, which is advantageous for trajectory planning and motion control of the moving platform.

3 Ay

 3 as a function of the coordinate O'(x,y,z) of the moving platform.

  are three constraint equations as below.

  In summary, when the coordinates of point O' on the moving platform 1 are known, each input values 1

FIGURE 6

 6 FIGURE 6 EXAMPLE OF PARALLEL SINGULARITY CONFIGURATION

  length of the link set by the PM, during the movement of the PM, there is always    , that is,

1 k and 2

 2

FIGURE 7 Fig. 9

 79 FIGURE 7 WORKSPACE AND THE PARALLEL SINGULARITY SITUATION

  FIGURE 9 DIFFERENT X-Y CROSS-SECTIONS IN THE WORKSPACE

ACKNOWLEDGMENTS:

This research is sponsored by the NSFC (Grant No.51475050 and No.51375062) and Jiangsu Key Development Project (No.BE2015043).

The direct solutions of the PM is calculated by MATLAB, as shown in Table 1. 18) and [START_REF] Zeng | Tri-pyramid Robot: Design and kinematic analysis of a 3-DOF translational parallel manipulator[END_REF], there are theoretically 32 inverse solutions, but for this example, 24 solutions belong to the virtual number solutions, which does not exist for the real robot. Therefore only 8 numerical solutions exist.

TABLE 2 THE VALUES OF INVERSE SOLUTIONS

No.


Xiaorong Zhu

Huiping Shen

Chengqi Wu

Damien Chablat
email: damien.chablat@cnrs.fr

Tingli Yang
email: yangtl@126.com

AUTOMATIC MOBILITY ANALYSIS OF PARALLEL MECHANISMS: AN ALGORITHM APPROACH BASED ON POSITION AND ORIENTATION CHARACTERISTIC EQUATIONS

published or not. The documents may come

INTRODUCTION

Mobility analysis is the basis for the research on mechanical synthesis, kinematics and dynamics analysis. However, it is often difficult due to the fact of dealing with linear dependency for complex parallel mechanisms (PMs). With the development of computer technology, it would be very helpful to develop software for the mobility analysis. To evaluate the PMs mobility, several methods have been developed based on the screw theory [START_REF] Huang Zhen | Advanced Spatial Mechanism[END_REF][START_REF] Kong | Type Synthesis of 3-DOF Parallel Manipulators Based on Screw Theory[END_REF], the Lie group [3~5], the Position and Orientation Characteristics (POC) theory [START_REF] Yang Tingli | Theory and Application of Robot Mechanism Topology[END_REF][START_REF] Shen Huiping | Position and Orientation Characteristic Based Method and Enlightenment for Topology Characteristic Analysis of Typical Parallel Mechanisms and Its Application[END_REF] and the geometric algebra [START_REF] Xinxue | Research on the Aanalysis Method of DOF of Parallel Mechanism in the Framework of Geometric Algebra[END_REF].

The approach upon screw theory has been successfully applied in many PMs. The calculation mainly involves the linear solution of screws which has the potential advantage of automatic analysis. However, it is difficult to obtain screws automatically. By establishing a coordinate system for a reference leg, Cao Wen'ao [START_REF] Cao Wenao | Digital Type Synthesis Theory of Spatial Multi-Loop Coupled Mechanisms[END_REF] realized an automatic analysis of the mobility of PMs.

The mobility analysis of the Lie group is based on the multiplication and intersection of displacement subgroup/submanifold. There are too many rules involved (over 107 rules) [START_REF] Rico | Lie Algebra and the Mobility of Kinematic Chains[J][END_REF], and not suitable for programming. The other related efforts can also be found in recent works [START_REF] Liao Ming | Symbolic Derivation of Position and Orientation Characteristics of Parallel Mechanisms[END_REF][START_REF] Ding | Computer-aided Structural Synthesis of 5-DOF Parallel Mechanisms and the Establishment of Kinematic Sstructure Databases[END_REF]. However, most of the existing methods rely on manual and are inefficient. So, the study of automatic mobility analysis can provide an effective and practical means for designers.

The mobility analysis method based on the POC theory has a clear formula and judgment criteria that are easy to use and program [START_REF] Yang Tingli | Theory and Application of Robot Mechanism Topology[END_REF]. To solve these problems, based on POC theory, a general computer-assisted platform for the analysis of PMs mobility will be written in this document. A digital model is proposed for topological structure that allows the POC to be mapped. An algorithm is established for calculation of POC of legs and PMs. Then, the principle and algorithm of PM mobility analysis are studied and software is developed. Finally, typical examples of parallel robots are used to show the efficiency of the software.

THEORY OF POC-BASED

Definition and related equations of POC

Definition of POC

To describe the relative motion characteristics between any two components in a mechanism, the POC is defined in [START_REF] Yang Tingli | Theory and Application of Robot Mechanism Topology[END_REF] as t1 r1 t (dir.) M= r (dir.)

     (1)
ξ = ξt + ξr ≤ F (2) where M is the POC, ξ is rank of the POC (i.e. number of independent elements), dir means the direction of output axis, F is the DOF of the mechanism.

POC equation for serial mechanisms

The POC of a serial mechanism is

m Li Ji i=1 M = M U (3)
where M Ji is the POC of the i th kinematic joint, m is the number of kinematic joints, ∪ is "Union" operation [START_REF] Yang Tingli | Theory and Application of Robot Mechanism Topology[END_REF].

POC equation for PMs

The POC of a PM is

1 pa Lj j1 MM     I (4)
where M pa is the POC of the moving platform, ν is the number of independent loops, M Lj is the POC of the j th leg for the same base point O', ∩ is "intersection" operation [START_REF] Yang Tingli | Theory and Application of Robot Mechanism Topology[END_REF].

Equation of mobility

A PM with (v+1) legs can be considered as a combination of v independent loops (SLC). The structure composition of PMs based on SLC [START_REF] Yang Tingli | Theory and Application of Robot Mechanism Topology[END_REF] is shown as Fig. 1. Two legs are chosen to form the first independent loop (SLC 1 ), and the mobile platform is considered an active component of the loop; Then, SLC 1 is regard as a whole (an equivalent sub-PM), and combined with another leg to form the second independent loop (SLC 2 ); and in the same way, (v+1) th leg combined with SLC (v- 1) to form the v th independent loop (SLC ν ).

FIGURE 1. THE BASIC INDEPENDENT LOOPS

COMPOSED OF PMS From this point, the number of mobility, i.e., the degree of freedom(DOF) of a PM can be calculated by Eq.( 5), and the POC of the moving platform can be used to convey the property of mobility.

j j m 1 i L j L i 1 j 1 j=1 j 1 F J = f                (5)
where,

j i ( j 1) j L b b i1 dim.(( M ) M )    I U (6)
ν = m-n+1 (7) where F is the DOF of a PM, J i is the DOF of the i th joint, f j is the DOF of all joints in the j th leg, n is the number of components, ξ Lj is the number of independent displacement equations of the j th loop,

Key technologies for automatic mobility analysis

① The topological description of PMs is one of the most significant issues. It should include complete topological information of PMs with the most concise representation so that it can be recognized and extracted automatically.

② POC description for a leg and for a PM is a fundamental problem in automatic mobility analysis. Although the symbolic description as in Eq.( 1) is good for expressing the geometric meaning, it is difficult to realize the transitivity of joint orientation relation.

③

The most critical issue is how to establish an algorithm for POC. In essential, the calculation of POC is to deal with linear dependency in complex PMs.

REPRESENTATION OF MECHANISM AND ITS POC

Description of the legs and their POC

Only simple legs constituted by single-DOF joints R and P are considered here. In general, joints in a PM are denoted by special symbols: R for revolute joint and P for prismatic joint.

Equivalent sub-PM

In a leg, geometric relationship between joints is divided into six types: parallel, vertical, coaxial, spatial common point, coplanar and arbitrary, denoted by "⊥ ", "||", "/", "*", "#" and "-" respectively.

Digital matrix of leg

A kind of ordered topology matrix (L) is developed for representing the topological structure of a leg in a PM. The diagonal element J i represents the type of pais from the fixed platform to the moving platform. Non-diagonal element N ij denotes the geometric relationship between joint i and joint j. The ordered topology matrix (L) is

1 1i 1f i1 i if f1 f i f J L J J NN N N NN                LL M O M M O M LL (8)
For convenience of programming, N ij (6 kinds of geometric relationship above-mentioned) is transformed into number 1~5, or 0 respectively. J i (R, P) is represented by 8 or 9. Thus, a leg can be expressed as a decimal matrix. Table 1 gives four typical legs and their matrices.

         La 0 0 1 M 1 1 0        (b) RRC leg (RRRP) 8 1 1 1 1 8 1 1 L 1 1 8 1 1 1 1 9            Lb 3 0 0 0 M 1 0 0 0        (c) PRRRR leg 9 2 2 1 1 2 8 1 2 2 L 2 1 8 2 2 1 2 2 8 1 1 2 2 1 8                Lc 3 0 0 0 0 M 0 1 0 1 0        (d)UPS leg (RRPRRR)
8 2 0 0 0 0 2 8 0 0 0 0 0 0 9 0 0 0 L 0 0 0 8 2 2 0 0 0 2 8 2 0 0 0 2 2 8

                 Ld 3 0 0 0 0 0 M 3 0 0 0 0 0       

Definition of POC matrix of a leg

The POC matrix of a leg should include not only the dimension of output, but also the direction of output. So, a kind of 2×f matrix (M) is proposed to map the POC of the leg to its topology matrix, that is

1 f 1 f t t M r r        L L (9)
The rules of the POC matrix of leg are as follows ① f is the number of the single-DOF joints in a leg.

② t i or r i (i=1~f) denotes translation and rotation output respectively, and values is 0, 1, 2 or 3. ③ t i or r i ≠0 indicates the existence of independent motion. ④ when t i or r i =3, the dimension of translation/rotation output is 3, and the direction is arbitrary in space and need not to be specified.

⑤ The dimension of translation/rotation is n ri i1 =r    / n ti i1 =t  
 , and the total dimension of output at the end of the leg is ξ=ξ r +ξ t . ⑥ Output direction of column i is related to the axes of the i th joint in the leg. a. t i =1 means there is an independent translation along the axis of the i th joint (r i =0) or along the normal plane of the i th joint (r i =1). b. when t i is 2, there are two independent translation in the normal plane of the i th joint axis.

POC matrix of sub-chain

Planar sub-chain (2-DOF G 2 and 3-DOF G 3 ) or spherical sub-chain (2-DOF S 2 and 3-DOF S 3 ) constitute a leg or part of a leg. Based on the definition of POC matrix, there are 7 types matrices for 15 sub-chains, as shown in Table 2.

TABLE 2 15 TYPES OF SUB-CHAINS AND THEIR POC No. Type Symbol Dimension Topology

Matrix POC matrix 1 G 2 R||R 2 8 1 1 8       1 0 1 0       2 R⊥ P 2 8 2 2 9       1 0 1 0       3 P⊥ R 2 9 2 2 8       01 01       4 G 3 R//R//R 3 8 1 1 1 8 1 1 1 8         200 1 0 0       5 R//R⊥ P 3 8 1 2 1 8 2 2 2 9         200 1 0 0       6 P⊥ R//R 3 9 2 2 2 8 1 2 1 8         0 2 0 0 1 0       7 R(⊥ P)//R 3 8 2 1 2 9 2 1 2 8         200 1 0 0       8 R(⊥ P)⊥P 3 8 2 2 2 9 2 2 2 9         200 1 0 0       9 P(⊥P)⊥ R 3 9 2 2 2 9 2 2 2 8         0 0 2 0 0 1       10 P(⊥ R)⊥P 3 9 2 2 2 8 2 2 2 9         0 2 0 0 1 0       11 S 2 R-R 2 8 0 08       0 0 1 1       12 R⊥ R(U) 2 8 2 28       0 0 1 1       13 S 3 R 1 -R 2 -R 3 3 8 0 0 0 8 0 0 0 8         000 111       4
No. Type Symbol Dimension Topology Matrix POC matrix 14 R 1 ,R 2 ,R 3 (S) 3 8 5 5 5 8 5 5 5 8

        000 111       15 R 1 ⊥ R 2 ⊥ R 3 ;R 1 ⊥ R 3 3 8 2 2 2 8 2 2 2 8         000 111      
POC Matrix of legs Take the leg UP as an example, it contains R 1 ⊥ R 2 ⊥ P 3 . Dimension of the output is 3 (2-rotation and 1-translation), and the rotation is around R 1 axis and R 2 axis, and the translation along P 3 axis. Thus, POC matrix of this leg is La

0 0 1 M 1 1 0       
.

Description for PMs and its POC

Description for PMs To make information concise, a new structural representation of PMs is proposed, which consists of legs and joints on the two platforms.

Firstly, legs of a PM are labeled in turn with the numbers 1~6. Then, all legs are represented using the ordered topology matrix L. Finally, the joints on two platforms are regarded as two correspond virtual legs, expressed by the ordered topology matrices M and B. Therefore, the PM can be described as PM=(L1, L2, ..., L(v+1), M, B) (10) The PM in Fig. 2

                 , 4 8 2 2 L 2 8 2 2 2 9         
Matrix M and B are 8 0 0 0 0 8 0 0 M 0 0 8 0 0 0 0 8

           , 8 0 0 0 0 8 0 0 B 0 0 8 0 0 0 0 8           

POC matrix of PMs

According to the intersection property of sets, the direction of translation/rotation output of a PM is related to one or more joints in a certain leg. Thus the No. of this leg should be recorded. So, a column matrix is defined as

1 n 1 pa 1 n 2 t t l M r r l              L L , (n≤6) (11)
where, l 1 and l 2 are the No. of the leg specifying the direction of translation/rotating output, n is the dimension of independent output. The PM shown in Fig. 2(A) has 2-rotation and 1translation. The POC matrix is pa

0 0 1 4 M 1 1 0 4              , the double 4
indicate that the PM has a translation along P 43 joints, and the directions of rotation are around R 41 and R 42 , both in the fourth leg.

AUTOMATIC GENERATION OF POC AND DIMENSION OF INDEPENDENT DISPLACEMENT EQUATIONS

Automatic POC generation of legs

According to Eq.( 3), POC of a leg is the union operation of the POC of all joints. Each leg can be considered as a serial of independent sub-chains and single-DOF joints. Thus, the subchains are extracted, and union operation of POCs are mapped to logical "OR" operation of matrices.

For convenience of programming, the POC matrix with dimension less than 6 is supplemented to the same dimension of 2×6 Thus, the POC of a leg can be computed by

3 2 3 2 L G G S S J POC M POC M M M M M       ( ) ( ) (13)
where, M G3 , M G2 , M S3 , M S2 and M J are the supplemented forms of the POC matrices of G 3 , G 2 , S 3 , S 2 and single-DOF joints.

Here, based on the topology matrices of legs and the extracted sub-chains and single-DOF joints, the POC of the leg can be generated as follows ① Recognize topology matrix of the leg, and extract G 3 , G 2 , S 3 , S 2 successively. ② Write out the supplemented matrix of G 3 , G 2 , S 3 , S 2 and single-DOF joints. ③ Carry out "OR" operation on supplemented matrices of G 3 , G 2 , S 3 , S 2 and single-DOF joints, and obtain the initial POC matrix of the leg, which denoted as Table 1. gives the matrices of four typical legs. Here take the leg P⊥ R//R⊥ R//R as an example.

① Recognize and sequentially extract planar and spherical sub-chains: one G 3 sub-chain P⊥ R//R and one G 2 planar sub-chaines R//R. ② Generate the supplemented POC matrix of G 3 and G 2 as G3 0 2 0 0 0 0 M 0 1 0 0 0 0

       , G2 0 0 0 1 0 0 M 0 0 0 1 0 0       
.

③ Carry out "or" operation of M G3 and M G2 , and obtain the POC:

L 0 2 0 1 0 0 3 0 0 0 0 0 M = 0 1 0 1 0 0 0 1 0 1 0 0             
, ξ r 0 =2, and ξ t 0 =3. ④ ξ t 0 = 3 and ξ r 0 = 2 ,so M L = M L 0 , that means

L 3 0 0 0 0 0 M 0 1 0 1 0 0       
, and ξ t = 3 indicates that the leg has three independent translation outputs, ξ r =2 indicates that it has two rotation outputs, r 2 =r 4 =1 means the two rotations are around the axes of R 2 and R 4 respectively.

Automatic POC generation of PMs

According to Eq.( 4), the POC of a PM is the intersection operation on all legs in the PM. Among them, the translational output of the PM is the intersection of translational outputs of all legs, and the same goes for the rotating outputs.

Given two legs L 1 and L 2 , numbered as l 1 and l 2 , the POC matrices are For convenience of description, the translational output of leg i is denoted as G i =(ξ ti , e i ) and the rotating output as H i =(ξ ri , s i ), where e i /s i are the direction of translational/rotating output of this leg. If ξ ti =0/ξ ri =0, e i /s i are Φ(empty set); if ξ ti =1/ξ ri =1, e i /s i are spatial lines; if ξ ti =2/ξ ri =2, e i /s i are plane; if ξ ti =3/ξ ri =3, e i /s i are "-".

i1 i2 i3 i4 i5 i6 Li i1 i2 i3 i4 i5 i6 t t t t t t M r r r r r r        (i=1 ,
Suppose the intersection result is denoted as

1 2 3 4 5 6 L1 L2 1 2 3 4 5 6 t t t t t t M r r r r r r        I , the dimension of translational/ rotating output is 6 i i1 r   / 6 i i1 t   respectively.
The translational output is G (1∩2) =(ξ t , e) and the rotating output is H (1∩2) =(ξ r , s), where e and s are the direction correspond.

Note that the legs L 1 and L 2 , which participate in computing, can be either the legs constituting the PM or the sub-PM composed by several legs. The POC of PMs can be calculated as follows.

Algorithm for translational output of PMs

The intersection operation of translational output of two legs is to solve G (1∩2) =G 1 ∩G 2 , which can be obtained according the rules shown in Table 3.

TABLE 3 INTERSECTION RULES FOR TRANSLATION L 1 L 2 ξ t1 =0 ξ t1 =1 ξ t1 =2 ξ t1 =3 ξ t2 =0 G 2 G 2 G 2 G 2 ξ t2 =1 G
ξ t2 =3 G 1 G 1 G 1 G 1

Algorithms for rotating output of PMs

Similarly, the intersection of rotating output of two POC matrices is to solve H (1 ∩ 2) =H 1 ∩H 2 . The intersection rules are listed in Table 4.

TABLE 4 INTERSECTION RULES FOR ROTATION L 1 L 2 ξ r1 =0 ξ r1 =1 ξ r1 =2 ξ r1 =3 ξ r2 =0 H 2 H 2 H 2 H 2 ξ r2 =1
ξ r2 =3 H 1 H 1 H 1 H 1
Accordance to the analysis above, the flow of generating POC matrices of PMs are as follows ① Input the number of legs, and POC of each leg. ② Intersection operation on translational outputs. ③ Intersection operation on rotating outputs. ④ Output the result of POC of the PM.

Here take the PM of 3-RRC as an example. As shown in Fig. 2(B). The topology matrices of legs, joints on the two platforms are described previously. ① Input the number of legs, and POC of each leg.

Number of legs is 3, and topology matrices of legs are i Li i 2 0 0 0 0 0 0 0 0 1 0 0 3 0 0 0 0 0 G M = 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 H

                           (i=1,2,3).
Thus, ξ ti =3 means arbitrary translation in space, and ξ ri =1 is one-dimension rotation around the R i1 axis in i th leg, i.e. s i =R i1 . ② Intersection operation on translational outputs According to

pa =M pa(1∩2) ∩M b3 = 2 G    φ ∩ 3 3 G H    = 3 G    φ
The output of 3-RRC PM is 3-translation, and the result is consistent with those in Ref [START_REF] Herve | The Lie Group of Rigid Body Displacements, a Fundamental Ttool for Mechanism Design[END_REF].

Calculating the number of independent displacement equations

Eq [START_REF] Yang Tingli | Theory and Application of Robot Mechanism Topology[END_REF] shows that solving the number of independent displacement equations of an independent loop (SLC) involves the union operation of POC matrices. For the legs L 1 and L 2 mentioned above, supposing that the result matrix of union is

1 2 3 4 5 6 L1 L2 1 2 3 4 5 6 t t t t t t M r r r r r r        
The dimension of independent output is

ξ Lr = ξ Lt =
Translational output is G (1∪2) = (ξ Lt , e) and rotating output is H (1∪2) = (ξ Lr , s), where e/s is the direction of translation/rotation. Then the number of the independent displacement equations of the independent loop is ξ L = ξ Lt + ξ Lr .

Algorithms for the dimension of translation output of SLC

Solving the dimension of translational output of a SLC is essentially to calculate ξ Lt =dim(G 1 ∪ G 2 ), whose operation rules are shown in Table 5.

TABLE 5. ALGORITHMS FOR THE DIMENSION OF TRANSLATIONAL OUTPUT

L 1 L 2 ξ t1 =0 ξ t1 =1 ξ t1 =2 ξ t1 =3 ξ t2 =0 0 1 2 3 ξ t2 =1 1 e 1 ||e 2 ξ Lt =1 e 1 ||e 2 ξ Lt =2 3 e 1 ||e 2 ξ Lt =2 e 1 ||e 2 ξ Lt =3 ξ t2 =2 2 e 1 ||e 2 ξ Lt =2 e 1 ||e 2 ξ Lt =2 3 e 1 ||e 2 ξ Lt =3 e 1 ||e 2 ξ Lt =3 ξ t2 =3 3 3 3 3
Similarly, the dimension of rotating output of a SLC is to calculate ξ Lr =dim(H 1 ∪ H 2 ), and its operation rules are shown in Table 6.

TABLE 6. ALGORITHMS FOR THE DIMENSION OF ROTATING OUTPUT L 1 L 2 ξ t1 =0 ξ t1 =1 ξ t1 =2 ξ t1 =3 ξ t2 =0 0 1 2 3 ξ t2 =1
ξ t2 =3 3 3 3 3

PROCEDURE OF AUTOMATIC MOBILITY ANALYSIS OF PMS

Based on the algorithm for POC of legs and PMs and for the number of the independent displacement equations, the mobility analysis flow for PMs is established as shown in Fig. 3.

Step1. Input topological structure of the PM.

Input topology matrices L 1 , …, L (ν+1) and number leg in sequence automatically, and input the axis relation matrix M and B on the moving platform and fixed platform. Obtain the number of the joints f i and the number of ν+1 legs. Step2. Calculate (ν+1) POC matrix of legs: M L1 ,...M L(ν+1) .

Step3. Calculate the sum of all joints in the PM:

f (1~(ν+1)) = +1 i 1 f i    .
Step4. Calculate the number of the independent displacement equations of the 1 st loop SLC 1 : ξ 1 ① Calculate the dimension of the independent translational output of SLC 1 :

ξ L1t = dim(G 1 ∪ G 2 ). ② Calculate the dimension of te independent rotating output of SLC 1 : ξ L1r = dim(H 1 ∪ H 2 ). ③ Calculate ξ L1 = ξ L1t +ξ L1r .
Step5. Calculate the POC matrix MP (1∩2) of the sub-PM P [START_REF] Huang Zhen | Advanced Spatial Mechanism[END_REF][START_REF] Kong | Type Synthesis of 3-DOF Parallel Manipulators Based on Screw Theory[END_REF] composed by the 1 th and 2 th legs.

① Calculate the translational output G (1∩2) =G 1 ∩G 2 . ② Calculate the rotating output H (1∩2) =H 1 ∩H 2 . ③ Get 12 (1 2) 12 G MP H        () () .
Step6. Calculate the POC matrix M pa(1~j) = M pa(1~(j-1)) ∩M bj of the sub-PM P (1-j) composed by the front j th (j=3,...,ν) legs. Step7. Calculate the number of the independent displacement equations of the j th (j=2,...,ν) independent loop SLC j :

ξ j =dim(M pa(1~j) ∪ M b(j+1) ).
Step8. Calculate the number of DOF:

  ν i 1~(ν+1) i=1 F=f -ξ  .
Step9. Determine the property of mobility of this PM:

M pa =M pa(1-ν) ∩M L(ν+1) . 6 i i1 r   6 i i1 t   FIGURE 3 THE FLOW CHART OF MOBILITY ANALYSIS BASED ON POC

SOFTWARE IMPLEMENTATION AND CASE STUDY

Software implementation

According to the basin principle discussed above, an automatic mobility analysis platform for PMs was developed by VC++ 6.0. Using the platform, large number of PMs listed in [START_REF] Yang Tingli | Theory and Application of Robot Mechanism Topology[END_REF] have been analyzed automatically based on the input of topology matrix and results have proved the validity of the software.

The input information mainly includes the number of legs and topology matrix. Using decimal array and easy to operate. The analysis results include the DOF, the dimension of the independent translational/rotating output of the moving platform, and their corresponding axis direction.

Case study

Detailed examples of mobility analysis of two PMs are presented in this section.

Mobility analysis of the Tricept PM

The PM of Tricept shown in Fig. 2(A) has 3 DOF which is one translation and two rotation. The legs are 3-UPS and 1-UP. The axes of joints on the two platforms are of general.

(1) Input topological structural matrices of the PM Topology matrix of UPS leg is i 8 2 1 0 0 2 9 2 0 0 L 1 2 8 0 0 0 0 0 8 5 0 0 0 5 8

               , (i=1~3).
Topology matrix of UP leg is 4

8 2 2 L 2 8 2 2 2 9          .
To the two virtual legs on the platforms, topology matrices M and B are 8 0 0 0 0 8 0 0 M 0 0 8 0 0 0 0 8

           , 8 0 0 0 0 8 0 0 B 0 0 8 0 0 0 0 8           
.

Number of legs ν+1=4 (2) Calculate POC matrix of legs M L1 ,...M L4 .

① i Li i 3 0 0 0 0 0 G M = 3 0 0 0 0 0 H             
, ξ ir =3, and ξ it =3 (i=1~3)

② 4 L4 4 0 0 1 0 0 0 G M = 1 1 0 0 0 0 H             
,and ξ 4t =1, ξ 4r =2.

(3) Get the sum of all joints in this PM: f (1-4) =6+6+6+3=21.

(4) Calculate the number of the independent displacement equations of the first loop (SLC 1 ) ① ξ ti =3, as shown in Table 5, ξ L1t =3, G (1∪2) =G 2 (i=1,2).

② ξ ri =3, then ξ L1r =3, H (1∪2) =H 2 (i=1,2).

③ So, ξ L1 =ξ L1t +ξ L1r =3+3=6.

(5) Calculate the POC matrix MP (1∩2) of sub-PM constituted by 1 st and 2 nd legs. ① ξ ti =3, as shown in Table 3, G (1∩2) =G 2 (i=1,2).

② ξ ri =3, as shown in Table 4, H (1∩2) =H 2 (i=1,2).

③ So, Bit "or" operation of the supplement matrix;

G G G = H H H                   3 ∩ . ( 7
       v i Li v i i f 1 1 1 F  M pa =M pa(1~ν) ∩M b(ν+1)
Calculating the POC M bi (i=1,…,ν+1) of leg.

Similarly, the number of the independent displacement equations of SLC 3 is ξ L3 =ξ L3t +ξ L3r =3+3=6.

(8) Get the number of DOF:

  3 i 1~4 i1 F=f - =21-6+6+6 =3    () .
(9) Get the property of mobility:

M pa =M pa(1-3) ∩M b4 = 3 G H       3 ∩ 4 G H       4 = 4 G H       4 .
The results show that this mechanism has 1T2R output, the rotating directions are around R 41 and R 41 , and the translational direction along P 43 joint. The result generated by software is shown in Fig. 4, which is consistent with those in Ref. [START_REF] Yang Tingli | Theory and Application of Robot Mechanism Topology[END_REF].

i 8 1 1 1 1 8 1 1 L 1 1 8 1 1 1 1 9            , (i=1~3).
To two virtual legs on the platforms, topology matrices M and B are

FIGURE 5. MOBILITY ANALYSIS OF 3-RRC PM

CONCLUSIONS

A matrix description for mapping the topological structure with POC of legs and PMs was established. It includes not only the type and size of the translational/rotary output, but also represents the orientation of the output axis. By extracting the planar and spherical sub-chains orderly, the POC of a leg can be transformed into the logical "OR" operation of the matrices. An algorithm for POC of a PM was established without manual intervention. Algorithm of mobility analysis automatically of PMs is proposed. Software for mobility analysis of PMs was created and typical examples were provided in detail to show its effectiveness.

  of the equivalent sub-PM composed by the front j legs, M L(j+1) is the POC of the (j+1) th leg.

  FIGURE 2 TWO KINDS OF PMSThe topology matrix L i (i=1~3) for 3-UPS legs and matrix L 4 for UP leg are respectively

④

  Process the translational output generated by parallel revolute joints. a. If ξ r 0 >3, then ξ r =3 and ξ t =ξ t 0 +ξ r 0 -3, and the translation newly generated lies in the normal plane of the single-DOF or S 2 rotational joints. b. If ξ r 0 ≤ 3 and ξ t 0 ≤ 3, the final POC matrix is M L =M L 0 .

6 )

 6 Calculate the POC matrix of sub-PM composed by the front three legs: M pa(1~3) = M pa(1~2) ∩M b3

FIGURE 4 .

 4 FIGURE 4. MOBILITY ANALYSIS OF TRICEPT Mobility analysis of a 3-RRC PM Fig.2(B) shows a 3-RRC PM. It has three identical RRC legs connecting the two platforms, labeled with the number 1~3. The topological structure of this PM is PM=(L 1 ,L 2 ,L 3 ,M,B). The fixed platform is a triangle, and axes of the joints on the two platforms are all coplanar. (1) Input topology matrices of the PM

③④,( 3 )

 3 Calculate POC matrix of legs M L1 ,...M L3 . ① Recognize and sequentially extract planar and spherical sub-chains: one G 3 sub-chain R//R//R and one single-DOF P joint. ② Generate the POC supplemented matrices of G 3 Carry out "OR" operation on POC matrices of all legs Li 2 ξ t =3 and ξ r =1, then M L =M La 0 and ξ it =3 indicate that the leg has three translational output along arbitrary direction, ξ ir =1 indicates that it has one rotating output, r i1 =1 means that the rotation is around the axes of R i1 , i.e. s i =R i1 . Calculate the sum of all joints in this PM: f (1-3) =4+4+4=12. (4) Calculate the number of the independent displacement equations of the first loop (SLC 1 ) ① ξ ti =3, as shown in

6 ) 7 )( 8 )( 9 )

 6789 Calculate the POC matrix of PM constituted by the front three legs: M pa(1~3) =M pa(1~2) ∩M L3 = Calculate the number of the independent displacement equations of SLC 2 : ξ L2 =dim(M pa(1~2) ∪ M b3 ) ① ξ t3 =3, as shown in Table 3, ξ L2t =3. ② ξ r(1-2) =0, and ξ r3 =1, as shown in Table 4, ξ L2r =1. ③ Thus, ξ L2 =ξ L2t +ξ L2r =3+1=4. Get the number of DOF: Get the property of mobility: M pa =M pa(1-3) = 3 G    φ . It shows that this PM has 3 DOF, three translation. The result automatically generated shown in Fig.5. The result is consistent with those in Ref[5].

TABLE 1 .

 1 TOPOLOGY MATRIX OF LEGS

	No.	Schematic diagram	Topology Matrix	POC matrix
	(a) UP leg (RRP)		8 2 2 L 2 8 2 2 2 9

Table 3 ,

 3 ξ ti = 3, then ξ L1t = 3, and G (1∩2) = G 2 .③ Intersection operation on rotating outputsAccording to Table4, ξ r1 =ξ r2 =ξ r3 =1, s 1 =R 11 , s 2 =R 21 , s 3 =R 31 and s 1 ||s 2 ||s 3 , then ξ r =0, s=Φ.

	④ Output the result of POC of this PM
	M

  ) Calculate the number of the independent displacement equations of SLC 2 : ξ L2 =dim(M pa(1~2) ∪ M b3 ). L2r =3. ③ then, ξ L2 =ξ L2t +ξ L2r =3+3=6.

	Calculate M pa(1~j) ==M pa(1~(j-1)) ∩M bj
	Calculate ξ Lj =dim{M pa(1~j) ∪M b(j+1) }
	j=j+1
	N
	j>ν+1
	Y
	Output:
	F and M pa
	End

① ξ t(1∩2) =3, and ξ t3 =3, then ξ L2t =3. ② ξ r(1∩2) =3, and ξ r3 =3, then ξ

Table 5 ,Table 6 Table 3 ,

 563 Table 4

	Start
	Input: number of legs ν+1;
	Topology matrix L i of i th leg;
	Topology matrix M and B on two platforms.
	Calculate the sum of DOF of joints in PM
	Extract planar sub-chains;
	Get the POC supplement matrix of
	sub-chains, and single-DOF joints.
	Calculate ξ L1 =dim{M L1 ∪M L2 }
	j=2

Table 5 ,

 5 ξ L1t =3 (i=1,2). ② ξ ri =1, and s 1 ||s 2 , then, ξ L1r =2 (i=1,2). ③ So, ξ L1 =ξ L1t +ξ L1r =3+2=5.(5) Calculate the POC matrix MP (1∩2) of sub-PM constituted by the 1 st and 2 nd legs ① ξ ti =3, as shown in Table3, G (1∩2) =G 2 (i=1,2). ② ξ ri =1, and s 1 ||s 2 , as shown in Table4, H (1∩2) =Φ (i=1,2).

ACKNOWLEDGMENTS:

This research is sponsored by the NSFC (Grant No. 51475050).


Fathia Chabi

Stéphane Champmartin

Christophe Sarraf

Ricardo Noguera

B Maurel

Comparison of various hemodynamic models for applications to cfd in stented arteries


Stanislav Aranovskiy
email: stanislav.aranovskiy@centralesupelec.fr

Rosane Ushirobira
email: rosane.ushirobira@inria.fr

Denis Efimov
email: denis.efimov@inria.fr

Géry Casiez
email: gery.casiez@univ-lille1.fr

An adaptive FIR filter for trajectory prediction and latency reduction in direct Human-Computer interactions

Keywords: human-computer interaction, latency reduction, adaptive FIR filter, human motion prediction, touchscreen

come

An adaptive FIR filter for trajectory prediction and latency reduction in direct Human-Computer interactions

Introduction

Typically, human-computer interactions can be divided into two categories: direct and indirect interactions. In indirect interactions, the input device (e.g. a mice or a trackball) and the output device (e.g. a display) are separated. In contrast, for direct interactions, the input and the output devices are coupled together, and the input and the system output (observed by the user) share the same screen. Examples of direct interactions include smartphones, pads and touch screens.

For any type of interaction some time is required for the input device to process the user input, to transfer it to the operating system, then to the specific application, to the graphical layer, and finally to display the output. This route gives rise to the end-to-end latency, which is the time lag between the input action performed by the user and the reaction output displayed by the hardware. The diagram of this route in direct interactions is given in Fig. 1. For modern touch-screen devices, it is reasonable to expect the end-to-end latency about 60 ms or more, as measured by Ng et al. (2012a). As shown by [START_REF] Casiez | Characterizing latency in touch and button-equipped interactive systems[END_REF], even for high-end touch interfaces, end-to-end latency is at least around 50 ms in an application simply changing the background color, where most of the latency comes from the display. End-to-end latency can be reduced with higher input and output frequencies (see Table 6 by [START_REF] Casiez | Characterizing latency in touch and button-equipped interactive systems[END_REF]), but this comes at the price of more expansive hardware and higher power consumption.

Measuring the end-to-end latency and its variations is not a trivial task itself, and several studies were taken on this issue recently [START_REF] Bérard | Two touch system latency estimators: high accuracy and low overhead[END_REF]; [START_REF] Casiez | Looking through the Eye of the Mouse : A Simple Method for Measuring End-toend Latency using an Optical Mouse[END_REF]). However, the detrimental impact of latency to user performance has been known for a long time, see [START_REF] Mackenzie | Lag as a determinant of human performance in interactive systems[END_REF]. Direct interactions are more sensitive to latency, and [START_REF] Jota | How fast is fast enough?: a study of the effects of latency in direct-touch pointing tasks[END_REF] found that latency greater than 25 ms can significantly affect the user performance in touch dragging tasks. At the same time, Ng et al. (2012b) show that latency as small as 10 ms can be perceived in direct interactions.

There are two ways to reduce the impact of latency. The first is at the hardware level, e.g. by using more reactive and advanced elements and by making the signal flow as fast as possible. This approach has two drawbacks, the high cost of advanced components and the increased energy consumption, which play a significant role for portable devices.

The second way to reduce latency is at the software level by using latency reduction algorithms. From the control point of view, this problem can be formulated as a trajectory prediction or forecasting problem, and convenient prediction methods can be used. However, methods based on underlying dynamic models cannot be easily applied for latency reduction since dynamic models of user behavior are not typically available. There exist several dynamic models of human movements, e.g. the open-loop minimumjerk optimal model used by [START_REF] Lank | Endpoint prediction using motion kinematics[END_REF], and the recently proposed dynamic models for indirect human-computer interaction [START_REF] Aranovskiy | Modeling pointing tasks in mouse-based human-computer interactions[END_REF][START_REF] Varnell | Stability and robustness analysis for human pointing motions with acceleration under feedback delays[END_REF]). However, [START_REF] Aranovskiy | Modeling pointing tasks in mouse-based human-computer interactions[END_REF][START_REF] Varnell | Stability and robustness analysis for human pointing motions with acceleration under feedback delays[END_REF] study the problem of modeling of human movements for a specific type of human-computer interaction, namely pointing motions in indirect interaction, and do not propose any solution for movement forecasting. Moreover, no systematic parameter identification procedure is known for these models restricting their applicability to the latency reduction problem. In contrast, this research is devoted to the problem of the user's trajectory prediction for general motions in direct human-computer interaction and it is not based on a underlying dynamic model.

The lack of models motivates the use of model-free prediction methods. For instance, a trajectory prediction using Kalman filter for a chain of integrators was proposed by [START_REF] Wu | On latency compensation and its effects on head-motion trajectories in virtual environments[END_REF], and a method based on the first-order Taylor series was used by [START_REF] Cattan | Reducing latency with a continuous prediction: Effects on users' performance in direct-touch target acquisitions[END_REF], where the velocity was estimated using the two most recent position measurements. In the recent paper by [START_REF] Ushirobira | A forecasting algorithm for latency compensation in indirect humancomputer interactions[END_REF] a forecasting algorithm based on the Taylor series expansion was proposed, where the derivatives were estimated using either algebraic [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF]), or homogeneous finitetime [START_REF] Perruquetti | Finite-time observers: application to secure communication[END_REF]) differentiators. The shortcoming is that, as it was stated by [START_REF] Ushirobira | A forecasting algorithm for latency compensation in indirect humancomputer interactions[END_REF], high non-linearity and interrelations between differentiator parameters make the tuning of these algorithms rather complicated. Also, some model-free approaches motivated by Kalman filter, curve fitting, and heuristic considerations can be found in the patents [START_REF] Moussavi | Methods and apparatus for incremental prediction of input device motion[END_REF], [START_REF] Qingkui | Curve fitting based touch trajectory smoothing method and system[END_REF], and [START_REF] Kim | Mobile terminal and touch coordinate predicting method thereof[END_REF], respectively.

To summarize, the problem of latency reduction can be translated to the problem of estimating the user's trajectory or, equivalently, to the problem of predicting further points on this trajectory. To this end, the state of the art is to apply a model-free estimation/prediction algorithm, which performs a prediction based on the past measurements. In this paper, a novel modelfree frequency-domain based approach for human movements forecasting is proposed. It can be shown that the ideal predictor is a non-causal linear time-invariant operator, which obviously cannot be implemented. However, assuming that human movements can be sufficiently well approached with a finite number of low-frequency components, the implementable forecasting algorithm is constructed as a causal low-frequency approximation of the ideal one. As it is shown in this paper, such an approximation can be either computed analytically with a priori knowledge about movement characteristics or obtained numerically as an optimization task that can be efficiently solved given a sufficiently rich dataset. Next, an adaptive modification of the design is proposed, which is intended to adapt to the changes of users and/or movement types. Experimental results illustrate the applicability of the proposed forecasting algorithm and its performance compared to other methods.

Preliminary results of this study were reported in [START_REF] Aranovskiy | Frequency domain forecasting approach for latency reduction in direct human-computer interaction[END_REF]. Extending that work, this manuscript contains i) the proof of the forecasting accuracy bound (15), see Proposition 1, ii) the detailed description of the experimental setup and iii) comparison to the methods currently available in patents [START_REF] Moussavi | Methods and apparatus for incremental prediction of input device motion[END_REF] and [START_REF] Qingkui | Curve fitting based touch trajectory smoothing method and system[END_REF].

The rest of the paper is organized as follows. First, the problem statement is given in Section 2. Next, the main forecasting algorithm is proposed in Section 3, and its adaptive version is considered in Section 4. Experimental results are described in Section 5, and concluding remarks are given in Section 6.

Problem statement

Consider a movement along a single axis. Motivated by physical reasoning, it is assumed that the trajectory x(t) is smooth and bounded. The raw measurement x(t k ) is the coordinate x measured at the time instance t k for an integer k ≥ 0. Given the latency value t L > 0, the goal is to construct the forecasting algorithm

x(t k + t L ) = F (x(t k ), t k ) ,
where x(t k ) := x(t k ), x(t k-1 ), . . . , x(t 0 ) , and x(t k + t L ) is the estimate of x(t k + t L ). In other words, given the history x(t k ) up to the time instant t k , the goal is to forecast the future value x(t k + t L ). The forecasting algorithm F(•) should be constructed in such a way that the forecasting error e x := x-x is small in the sense of a certain norm.

It is also assumed that the coordinate x is measured with the constant sampling interval T s and t k = kT s . Moreover, the latency value t L is the integer number L > 0 of sampling intervals, t L = LT s , and t k + t L = t k+L .

The forecasting algorithm

An ideal predictor approximation

The ideal predictor in discrete-time domain is the L steps ahead time shift, i.e. F(x, t k ) = x(t k+L ) = q L x(t k ), where q is the one step forward time shift operator. The ideal predictor has the following linear transfer function

F id (z) = z L ,
where z is a complex variable. In the frequency domain, one obtains

F id (iω) = e iωt L ,
where i is the imaginary unit and ω is the frequency [START_REF] Goodwin | Control system design[END_REF]). It follows

|F id (iω)| ≡ 1, arg F id (iω) = ωt L .
Obviously, such a transfer function is not causal and cannot be implemented. However, this concept of the ideal predictor motivates the construction of the forecasting algorithm as an approximation of the ideal one in the frequency domain. Due to the assumptions on the trajectory x(t), for any time instant t k and any window size T W , with

t k > T W > 0, one can represent x(t) for t ∈ [t k -T W , t k ] as x(t) = ∞ i=0 A i sin(ω i t) + B i cos(ω i t), (1)
where A i and B i are some constants and

ω i = 2πi T W .
Next assume that there exists a known value ω † , such that any reasonable trajectory generated by a user in human-computer interaction can be sufficiently well approximated with the first N † + 1 elements of the series expansion (1), where

N † is chosen such that ω N † ≤ ω † and ω N † +1 > ω † ; this approximation is further denoted as x † , x † (t) = B 0 + N † i=1 A i sin(ω i t) + B i cos(ω i t). (2)
To support this assumption, note that a user cannot provide infinite accelerations or decelerations to his hand (wrist). Consider a curve plotted by the user in the x -y plane without lifting up the stylus. Even if this curve is not smooth, the trajectory x(t) is necessary continuous and sufficiently smooth.

Then the forecast x † (t k + t L ) of this approximation can be used as the forecast of the coordinate x(t k + t L ). The rest of the expansion (1)

δ x (t) := x(t) -x † (t) = ∞ i=N † +1 A i sin(ω i t) + B i cos(ω i t)
is bounded and represents both negligible high-frequency trajectory deviations and possible measurement noise [START_REF] Goodwin | Control system design[END_REF]). Assuming that the coefficients of the approximation (2) are slow-varying during the latency interval t L , i.e. almost constant, the forecast x † of the approximation x † is given by

x † (t k + t L ) = B 0 + N † i=1 A i sin(ω i (t k + t L )) + B i cos(ω i (t k + t L )).
This forecast can be constructed by applying a linear operator F ap to the signal x † , where F ap satisfies

F ap (iω i ) = F id (iω i ), i ∈ {0, 1, . . . , N † }.
(3)

The operator F ap is thus an approximation of the ideal predictor F id : its frequency response coincides with the frequency response of F id for frequencies ω i ≤ ω † , and it does not necessary coincide for others. However, since the approximation x † is not known, i.e. it is not measured, the operator F ap is applied to the measured signal x, and the frequency-domain approximationbased forecasting algorithm takes the form

x(t k + t L ) = F ap [x(t)] = x † (t k + t L ) + F ap [δ x (t)], (4)
where F ap [(•)(t)] denotes the operator F ap applied to the signal (•)(t). The term F ap [δ x (t)] in the right-hand side of (4) represents the high-frequency distorting components and its impact should be minimized when possible.

Approximation using a finite impulse response filter

For any fixed ω † the operator F ap can be chosen as a causal stable linear time-invariant system of order N ≥ 2N † . One particular choice widely accepted in digital signal processing is a Finite Impulse Response (FIR) filter

F ap (z) = c 0 + c 1 z -1 + . . . + c N z -N = N i=0 c i z -i yielding the forecasting algorithm x(t k + t L ) = c 0 x(t k ) + c 1 x(t k-1 ) + . . . + c N x(t k-N ).
Given N † + 1 distinct frequencies ω i , i = 0, 1, . . . , N † , where

ω N † < π T s , (5)
there always exists a unique FIR filter F ap of order N = 2N † satisfying (3), which can be constructed by solving a system of the following N + 1 linear equations:

     1 N +1 c = 1, Re (Ω) c = Re f id , Im (Ω) c = Im f id , (6)
where c := c 0 c 1 . . . c N , Re(•) and Im(•) are the real and the imaginary parts of a complex value (•), respectively, 1 N +1 is the vector of 1's of dimension N + 1, f id is a N † -dimensional vector with complex entries, 2013)) that the system (6) has a unique solution if for all i ∈ {1, . . . , N † } it holds 0 < ω i T s < π. The latter is satisfied due to (5). An example of frequency responses of such approximating FIR filters with N = 2N † for T W = 2s, L = 50 and T s = 1ms is given in Fig. 2 for ω † = 3Hz and ω † = 7Hz.

f id := e iω 1 t L e iω 2 t L . . . e iω N † t L , 0
As it can be seen from Fig. 2, the filter computed in this way has significant amplification for high frequencies, and it increases impact of the F ap [δ x (t)] term, which can lead to unacceptable fast oscillations in the forecast x, known as jitter. One possible way to overcome this drawback is to over-parametrize the FIR filter, i.e. to choose N > 2N † . Obviously, in such a case there exists an infinite number of FIR filters satisfying the condition (3), and the predictor can be constructed via optimization procedure. For example, if no a priori information about δ x is known, then it is reasonable to choose c as c = arg min c F ap (z) ∞ subject to (3).

Otherwise, if some information about δ x is available, e.g. its frequency range, the cost function may be modified to take this information into account.

Data-based filter tuning

The tuning procedure proposed in Subsection 3.2 requires a priori knowledge of ω † , as well as δ x frequency range, if available. However, these parameters are not easy to obtain: they might depend on a specific input device and thus are hardware-dependent. That makes the analytical tuning questionable and gives rise to the experiment-based tuning approach.

Assume that a sufficiently large trajectory record R is available representing all kinds of possible/admissible movements. The record R consists of N R measurements

R := {x(t k ), k = 1, 2, . . . , N R }.
Given the FIR filter order N the forecasting error for k > L + N is computed as

e x (t k ) = x(t k ) -x(t k ) = c 0 x(t k-L ) + c 1 x(t k-L-1 ) + . . . + c N x(t k-L-N ) -x(t k ) = -x(t k ) + N i=0 c i x(t k-L-i ).
Define the vector of the forecasting errors as

ēx := e x (t N R ) e x (t N R -1 ) . . . e x (t L+N +1 ) .
Then the FIR filter coefficients c can be tuned in order to minimize a certain cost function J(ē x ). Particularly, if the cost function is chosen as a vector norm

J(ē x ) = ēx p , p = 1, 2, ∞,
then the tuning procedure is translated to the linear programming problem for p = 1, ∞ or to the least squares problem for p = 2, which can be efficiently solved even for a relatively large record size N R [START_REF] Sturm | Using sedumi 1.02, a matlab toolbox for optimization over symmetric cones[END_REF]; [START_REF] Currie | Opti: Lowering the barrier between open source optimizers and the industrial matlab user[END_REF]).

It is expected that the operator F ap tuned with a sufficiently large record R represents the best possible approximation of the ideal predictor F id with respect to the ω † and δ x information explicitly contained in the record.

Adaptive forecast

The approach considered in Subsections 3.2 and 3.3 proposes to design a single fixed-gain FIR filter and to use it for forecasting regardless the exact movement performed by a user. It is assumed that this fixed FIR filter is able to forecast sufficiently well for all possible trajectories and for all users. On the other hand, it is reasonable to assume that the ω † parameter and the characteristics of δ x (t) can vary for different users and/or for different kinds of movements, e.g. in drawing or pointing tasks. In this case, the FIR filter tuned for a specific user and/or movement can provide a better forecast accuracy for these specific conditions than a general fixed filter. In this section, on-line tuning of the FIR filter parameters is proposed in order to adapt to the exact user and/or to the specific movement being executed at the moment.

Recursive least-squares algorithm

To this end, rewrite the forecasting algorithm as

x(t k+L ) = φ (t k )c(t k ), (7)
where

φ(t k ) := x(t k ) x(t k-1 ) . . . x(t k-N ) (8)
and c(t k ) is the vector of FIR filter coefficients updated at the each step. The goal is to update the coefficients c(t k ) in order to minimize the cost function

J(t k ) = k i=L+N +1 λ k-i x(t i ) -φ (t i-L )c(t k ) 2 ,
where 0 < λ ≤ 1 is the forgetting factor, which provides higher weight to recent measurements. The desired update law is given with the recursive least-squares algorithm [START_REF] Aström | Adaptive Control, 2nd Edition[END_REF]):

g(t k ) = P(t k-1 )φ(t k-L ) λ + φ (t k-L )P(t k-1 )φ(t k-L ) , P(t k ) = 1 λ P(t k-1 ) -g(t k )φ (t k-L )P(t k-1 ) , c(t k ) = c(t k-1 ) + g(t k ) x(t k ) -φ (t k-L )c(t k-1 ) . ( 9
)
At each step k the new measurement x(t k ) is available and the gains c are updated using the earlier measurements φ(t k-L ). The updated gains are further used in the forecasting algorithm ( 7). The drawback of this approach is that the most recent measurements x(t k-1 ), . . . , x(t k-L+1 ) are not actually used in the gains update. In other words, if any changes in the movement have occurred, they do not impact the gains for the next L steps.

Adaptive predictor

To take into account the most recent measurements, the following procedure is proposed. First, the one-step-ahead adaptive predictor is constructed. Next, the gains of this predictor are used to compute the gains c. The onestep-ahead predictor is given by

x(t k+1 ) = φ (t k )θ(t k ), ( 10
)
where θ ∈ R N +1 is the gains vector. The gains θ can be adaptively tuned using the recursive least square algorithm ( 9) and substituting L = 1 and θ in the place of c.

Next, to compute c(t k ) it is assumed that θ(t k ) is slow-varying for a time window of at least L steps ahead, and the trajectory x can be approximated on this window with the following autoregressive model with time-varying gains

x(t k+1 ) = φ (t k )θ(t k ), ( 11
)
where for all i ∈ {0, 1, . . . , L} it holds

θ(t k+i ) -θ(t k ) ≤ θ , ( 12
)
and • is the Euclidean norm.

Proposition 1. Given the autoregressive model (11) under assumption (12), set the vector c(t k ) as

c(t k ) = θ(t k ) e 1 . . . e N L e 1 , (13)
where e i denotes the i-th Euclidean basis vector, i.e. the N + 1 dimensional vector with a 1 in the i-th coordinate and 0's elsewhere. Then it holds

x(t k+L ) = φ (t k )c(t k ) + σ x (t k ), (14)
and

|σ x (t k )| ≤ √ N + 1x ∞ H( θ , θ(t k ) ) θ , (15)
where the hardware-dependent value x ∞ is the maximum admissible position value, the function

H( θ , θ(t k ) ) := L-1 i=1 L -1 i i-1 θ ( θ(t k ) + 1) L-i , (16)
and n k is the binomial coefficient. Proof. Define the d-steps-ahead shift vector β(d, t k ) ∈ R N +1 such that for the autoregressive model ( 11) the following holds

x(t k+d ) = φ (t k )β(d, t k ). ( 17
) Obviously, β(1, t k ) = θ(t k ) and c(t k ) = β(L, t k ). Moreover, for all k > N it holds β(0, t k ) = e 1 , β(-1, t k ) = e 2 , and so on up to β(-N, t k ) = e N +1 . Next define a (N + 1) × (N + 1) matrix B(d, t k ) as B(d, t k ) := β(d, t k ) β(d -1, t k ) . . . β(d -N, t k ) ,
where by definition

B(1, t k ) = β(1, t k ) e 1 . . . e N . (18)
Recalling (8) it can be verified that

x(t k+d+1 ) = φ (t k+d )β(1, t k+d ) = x(t k+d ) x(t k+d-1 ) . . . x(t k+d-N ) β(1, t k+d ) = φ (t k ) β(d, t k ) . . . β(d -N, t k ) β(1, t k+d ) = φ (t k )B(d, t k )β(1, t k+d ).
On the other hand, from (17) it follows

x(t k+d+1 ) = φ (t k )β(d + 1, t k ),
and thus

β(d + 1, t k ) = B(d, t k )β(1, t k+d ). (19)
Note that by definition the first column of B(d+1, t k ) is β(d + 1, t k ), and the i-th column of B(d + 1, t k ) coincides with the (i -1)-th column of B(d, t k ) for i varying from 2 to N + 1. Then from (19) and recalling (18) one obtains

B(d + 1, t k ) = B(d, t k )B(1, t k+d ).
The latter implies

B(2, t k ) = B(1, t k )B(1, t k+1 ), B(3, t k ) = B(1, t k )B(1, t k+1 )B(1, t k+2 ), B(4, t k ) = B(1, t k )B(1, t k+1 )B(1, t k+2 )B(1, t k+3 ), and
B(d, t k ) = d-1 i=0 B(1, t k+i ). ( 20
) Therefore x(t k+d ) = φ (t k )B(d, t k )e 1 = φ (t k )β(d, t k ).
However, at the step k the only available value is θ(t k ), that gives β(1, t k ) and B(1, t k ), but not B(d, t k ). With the definitions in use, ( 13) can be rewritten as

c(t k ) = (B(1, t k )) L e 1 ,
and σ x (t k ) in ( 14) is thus

σ x (t k ) = x(t k+L ) -φ (t k ) (B(1, t k )) L e 1 = φ (t k ) B(L, t k ) -(B(1, t k )) L e 1 . Define ∆B(d) := B(1, t k+d ) -B(1, t k )
and note that

B(1, t k ) = sup x =1 B(1, t k )x = sup x =1 θ(t k ) e 1 . . . e N x = sup x =1 θ(t k )x 1 + x 2 x 3 . . . x N +1 0 ≤ θ(t k ) + 1,
and for all i ∈ {0, 1, . . . , L}

∆B(i) ≤ θ .
From (20) it follows

B(L, t k ) = B(1, t k ) L-1 i=1 (B(1, t k ) + ∆B(i)) and B(L, t k ) -(B(1, t k )) L = B(1, t k ) L-1 i=1 (B(1, t k ) + ∆B(i)) -(B(1, t k )) L-1 . Therefore B(L, t k ) -(B(1, t k )) L ≤ ( θ(t k ) + 1) L-1 i=1 L -1 i i θ ( θ(t k ) + 1) L-i-1 = θ H( θ , θ(t k ) ),
where H(•, •) is defined in ( 16). Noting that

φ(t k ) ≤ √ N + 1x ∞ ,
the upper bound (15) follows.

Remark 1. The upper bound (15) allows to establish some qualitative conclusions, e.g. a relation between the upper bound and the maximum admissible position x ∞ . One particularly important observation is that if θ(k) remains constant and θ equals zero, then |σ x (k)| is zero as well.

Finally, the adaptive forecasting algorithm that uses the most recent measurements can be formulated as follows.

Algorithm 1:

Step 1. At the k-th step get the new measurement x(t k ).

Step 2. Update the one-step-ahead predictor gains θ(t k ) using algorithm (9) with L = 1 and substituting θ in the place of c.

Step 3. Compute c(t k ) as ( 13) using the one-step-ahead predictor gains θ(t k ).

Step 4. Compute forecast x(t k+L ) as (7).

Experimental results

An experiment where three participants were asked to draw some predefined shapes on a screen using a stylus is considered. The shapes were divided into three groups: a) the geometric shapes given in Fig. 3,b) the Latin capital letters from A to F, and c) the digits from 0 to 9.

The hardware used in the experiment is the Touch X by 3D Systems (2019), previously known as Phantom Desktop, which is a pen-type haptic device. The Touch X is coupled with a horizontally oriented 15" display and calibrated to get the stylus tip position co-located with the display, using an optimal affine transformation between two 3D point sets using the RANdom SAmple Consensus (RANSAC) algorithm [START_REF] Fischler | Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography[END_REF]). Thus, this experimental setup mimics a touchscreen with a pen.

The measurements were captured with the sampling interval T s = 1ms, and the expected latency value t L = 50 ms, that is L = 50. The movements were performed in the plane, and the x and y axes are considered separately (human movements along different axes typically have different magnitudes and velocities, and the predictors for x and y may have different order N and operate over time windows of different length L).

This section presents a comparison of online forecasting algorithms which do not require a rich dataset to tune them. These algorithms are motivated by the patents [START_REF] Moussavi | Methods and apparatus for incremental prediction of input device motion[END_REF] and [START_REF] Qingkui | Curve fitting based touch trajectory smoothing method and system[END_REF] and are based on Kalman filtering and curve fitting. The 3-rd order Kalman filter is chosen, where the process and measurement model matrices are given by [START_REF] Moussavi | Methods and apparatus for incremental prediction of input device motion[END_REF]. For the curve fitting forecast, the 2-nd order polynomial fitted with the last 30 samples is used. These two predictors are compared with the adaptive FIR filter given by Algorithm 1. The FIR filters is of order N = 15, that implies 16 tunable gains. The forgetting factor is a tuning parameter and is chosen such that it provides the discounting factor 1 2 for the measurements obtained one second ago yielding where v(t k ) is the velocity estimate, and the differentiation window length d = 5. It was found that the majority of the movements is performed with the absolute velocity in the range form zero to 1200 pixels per second. This range is uniformly divided into 20 intervals, and for each interval the mean absolute error (MAE) is computed for the adaptive FIR, Kalman filter, and curve fitting predictors. The resulting plot is given in Figure 4. As it can be seen, the adaptive predictor outperforms it competitors in the considered velocity range.

λ 1 Ts = 2 ⇔ λ = exp(-T s log(2)) ≈ 0.9993.
Remark 2. As it is discussed in Section 3.3, the fixed-gains FIR filter requires sufficiently rich records to be tuned off-line, and thus it is not included in comparison. However, the frequency response of the fixed-gains FIR filter optimized over the whole dataset is given as an illustration in Fig. 5. The phase response shows that the ω † value approximated from the given experi- mental data is close to 5 Hz.

Conclusion

The problem of latency compensation in direct human-computer interactions was considered in this paper. This problem can be formulated as a forecasting algorithm design, and a novel approach was proposed, based on a frequency-domain approximation of the ideal predictor. Such an approximation can be either computed analytically or constructed by solving optimization problem for a sufficiently rich dataset. Also, an adaptive modification of the forecasting algorithm was developed allowing possible variations of the user behavior to be taken into account. Experimental studies illustrate the applicability of the proposed solution for a trajectory forecasting with reasonable accuracy.

One possible direction of further studies is to apply for adaptive parameters tuning some techniques with enhanced transient performance, e.g. the one recently proposed by Aranovskiy et al. (2017a). Hopefully, this modification can reduce the maximum prediction error value M E. Another challenging problem is to consider model-based trajectory prediction algorithms that involve the design of user behavior models.

Acknowledgments

This work was supported by ANR (TurboTouch, ANR-14-CE24-0009). Stanislav Aranovskiy is supported by the Russian Science Foundation grant

Conflict of interest

None declared.

Figure 1 :

 1 Figure 1: Direct interaction diagram (Ushirobira et al. (2016)).

  Magnitude (in dB) versus frequency (in Hz).

  Phase (in deg) versus frequency (in Hz).

Figure 2 :

 2 Figure2: Comparison of the ideal predictor and its approximations given by FIR filters for ω † = 3Hz and ω † = 7Hz

Figure 3 :

 3 Figure 3: Geometric shapes used in the experiment.

Figure 4 :

 4 Figure 4: MAE (in pixels) versus the movement velocity (in pixels per second).

  Phase (in deg) versus frequency (in Hz).

Figure 5 :

 5 Figure5: Comparison of the ideal predictor and its approximations given by fixed-gains FIR filter tuned with N = 15 using the whole dataset

Table 1 :

 1 Prediction performance comparison divided by figure groups.

		M AE M SE M E
		pixel pixel 2 pixel
		Geometric shapes
	Kalman filter 2.05	9.89 56.94
	Curve fitting	2.27	11.44 50.09
	Algorithm 1	1.97	8.22 46.82
		Latin letters
	Kalman filter 2.84	16.22 36.12
	Curve fitting	3.21	20.09 37.21
	Algorithm 1	2.33	10.48 39.76
			Digits
	Kalman filter 3.13	18.56 34.97
	Curve fitting	3.63	24.71 34.12
	Algorithm 1	2.32	10.13 25.77

Table 2 :

 2 Prediction performance comparison for the whole dataset.

		M AE M SE M E
		pixel pixel 2 pixel
	Kalman filter 2.70	15.15 56.94
	Curve fitting	3.08	19.17 50.09
	Algorithm 1	2.21	9.66 46.82

The performance indices used for comparison are:

• the mean absolute error M AE computed as

where N samples is the number of samples used for performance evaluation;

• the mean squared error M SE computed as

x (t k );

• the maximum error M E computed as

Note that MAE and ME are measured in pixels, while MSE is measured in pixels 2 .

Comparison results divided by figure groups are given in Table 1, and the same performance indices computed for the whole dataset are given in Table 2. It can be seen from the experimental results that the proposed adaptive FIR prediction algorithm provides the most accurate prediction in the terms of M AE, M SE and M E. The only case when the adaptive FIR has the M E value slightly higher than the competitors is the Latin Letters scenario and probably this can be related to transients in adaptive gains tuning. The Kalman filter predictor has the second best accuracy in terms of M AE and M SE. The M E value of the Kalman filter predictor is, in general, higher than the one of the curve fitting predictor, which, however, has the poorest M AE and M SE accuracy.

To evaluate how does prediction error depend on the movement velocity, the numerical differentiation was performed as v(t k ) =

x(t k+d ) -x(t k-d ) 2dT s ,