ABCM Symposium Series in Mechatronics - Vol. 3 - pp.316-325 c 2008 by ABCM Copyright °
A COLLISION AVOIDANCE METHOD USING ASSUR VIRTUAL CHAINS Henrique Simas,
[email protected] Universidade do Vale do Itajaí - Centro de Ensino São José Curso de Engenharia de Computação Rodovia SC 407, km 4 - Sertão do Imaruim 88122 000 - São José, SC. Brasil
Daniel Fontan Maia da Cruz,
[email protected] Raul Guenther,
[email protected] Daniel Martins,
[email protected] Universidade Federal de Santa Catarina Departamento de Engenharia Mecânica Laboratório de Robótica Campus Universitário - Trindade 88040-900 – Florianópolis, SC. Brasil
Abstract. In hydroelectric power plants, the rotor blades are eroded by the cavitation process. The erosion results in craters that are usually recovered by a manual welding process. The ROBOTURB project developed an automatized system, where a robot is used for recovery the rotor blade surfaces by welding process. The robot workspace is constrained by freeform surfaces and the collision imminency is constant. For this reason, the robot is redundant, making possible the collisions avoidance, keeping the endeffector trajectory tracking. In recent years, the collision avoidance problem has been dealt with algorithms based on artificial intelligence. This article considers a methodology for implementation a deterministic collision avoidance algorithm. The Solution of this proposal is based on the method of kinematic constraint and on the use of Assur virtual chains. The method consists in the obstacle identification and in the definition of a free collison workspace area. Thus, it is used a Assur virtual chain to detect and avoid the collision by evaluating the distance between a point at the robot and the obstacle. The Assur virtual chain is responsible for maintaining the robot inside the free collision workspace. This process is carried out for each point of the desired trajectory. Keywords: Differential kinematics, Assur virtual chains, Screw theory, Collision avoidance, Redundant robots 1.
INTRODUCTION
For many reasons, such as safety and economics, it is desirable to develop robots able to perform more complex tasks with minimal or no human intervention. A kind of these complex task is that in which the robot should avoid collisions autonomously. An example of this practical situation occurs in the Roboturb project, which aims to develop an automatized system for the repair of regions eroded by cavitation in the rotors of hydraulic turbines. The use of this system intends to eliminate the unhealthy work and to improve the welding quality in the repair. The central element of the system is a robot to operate in the confined space of the canal of the rotor, measuring the eroded region and performing the recovery welding. To reach this goal, the robot movements should be precise and the manipulator links cannot collide with the surfaces of the blades, i.e., the robot workspace is constrained. To make it possible to avoid collisions while the tasks are executed we designed a 7-DOF redundant robot (Guenther et al., 2000). This work was initially developed to solve the collision avoiding related to that project. Numerous methods have been proposed to achieve collision-free path planning, where the mainly application is for mobile robots. However, in the case of manipulator arms, the structure of the robot adds to the complexity of the solution (Soucy and Payeur, 2005). The numerous approaches for collision avoidance available on the literature are based basically on the potential fields' method, fuzzy logic or neural networks. Potential fields method is based on defining an attractive field, which directs the robot towards its objective, and a repulsive field that pushes the robot away from obstacles. Many researches have been developed using this approach (see (Muller, 2004), (Caselli et al, 2002) and (Piaggio and Sgorbissa, 1999)). However, in general these approaches do not offer a solution that can be used with generic robot architectures (Soucy, 2005). The methods based on fuzzy logic and neural networks are more complex methods, which are normally combined with potential fields concepts, as can be seen in ((Soucy and Payeur, 2005), (Soucy, 2005) and (Tian and Mao, 2002)). These approaches are considered as non-deterministic methodologies where it is not possible to proof the stability. The main contribution of this paper is to introduce a new method to avoid collision based on the use of kinematic constraints (Simas, 2005) and Assur virtual chains (Campos et al., 2005). To present the method we shortly reproduce three fundamental kinematic tools used in the paper: the theory of differential kinematics based on the screw representation; the Davies Method; and the Assur virtual chain concept (see (Davies, 1981) and (Bonilla, 2004) for more details). In the sequence, we present the proposed method and an application example based on Roboturb project robot.
2. FUNDAMENTAL KINEMATIC TOOLS Our approach is based on the screw representation of differential kinematics, on the Davies method, and on the Assur virtual chain concept, which are shortly presented in this section. 2.1. Screw representation of differential kinematics The Mozzi theorem states that the general spatial differential motion of a rigid body consists of a differential rotation about, and a differential rotation along a axis named instantaneous screw axis (see Ceccarelli (2000)). In this way the velocities of the points of a rigid body with respect to an inertial reference frame O-xyz may be represented by a differential rotation ω about the instantaneous screw axis and a simultaneously differential translation τ about this axis. The complete movement of the rigid body, combining rotation and translation, is called screw movement or twist and is here denoted by $ . Figure 1 shows a body “twisting” around the instantaneous screw axis. The ratio of the linear velocity and the angular velocity is called pitch of the screw h = τ ω .
Figure 1. Screw movement or twist. The twist may be expressed by a pair of vectors, i.e. $ = [ω T ; V pT ] , where ω represents the angular velocity of the T
body with respect to the inertial frame, and V p represents the linear velocity of a point P attached to the body which is instantaneously coincident with the origin O of the reference frame. A twist may be decomposed into its amplitude and its corresponding normalized screw. The twist amplitude, denoted as q& in this work, is either the magnitude of the angular velocity of the body, ω , if the kinematic pair is rotative or helical, or the magnitude of the linear velocity, V p , if the kinematic pair is prismatic. The normalized screw, $ˆ , is a twist in which the magnitude is factored out, i.e.
$ = $ˆq&
(1)
The twist given in Eq. (1) expresses the general spatial differential motion movement (velocity) of a rigid body with respect to an inertial reference frame O-xyz. The twist could also represent the movement between two adjacent links of a kinematic chain. In this case, the twist $i represents the movement of link i with respect to link (i − 1) . 2.2. Davies method Davies method is a systematic way to relate the joint velocities in closed kinematic chains. Davies (Davies, 1981, 2000) derives a solution to the differential kinematics of closed kinematic chains from the Kirchhoff circulation law for electrical circuits. The resulting Kirchhoff-Davies circulation law states that “The algebraic sum of relative velocities of kinematic pairs along any closed kinematic chain is zero” (Davies, 1981). We use this law to obtain the relationship among the velocities of a closed kinematic chain as in Campos et al. (2005) and in Santos et al. (2006). In this way, considering that the velocity of a link with respect to itself is null, the circulation law could be expressed as n
∑ $i = 0 i =1
(2)
where 0 is a vector which dimension corresponds to the dimension of the twist $i . According to the above introduced normalized screw definition this equation may be rewritten as n
∑ $ˆi q&i = 0
(3)
i =1
where $ˆi and q&i represent the normalized screw and the magnitude of twist $i , respectively. Equation (3) is the constraint equation which, in general could be written as Nq& = 0
(4)
[
]
where N = $ˆ1 $ˆ2 ... $ˆn is the network matrix containing the normalized screws which signs depend on the screw definition in the circuit orientation, and q& = [q&1 q&2 ... q&n ] is the magnitude vector. In case the kinematic chain has more than one closed kinematic chain, say k closed kinematic chains, the KirchhoffDavies circulation law is applied to each one of the chain closed circuits. In this case the network matrix has k rows of screws, each one corresponding to a closed circuit. A closed kinematic chain as actuated joints, here named primary joints, and passive joints, here named secondary joints. The constraint equation, Eq. (4), allows calculate the secondary joint velocities as functions of the primary joint velocities. To this end the constraint equation is rearranged highlighting the primary and secondary joint velocities. So, Eq. (4) can be written as follows:
[N
p
⎡q& p ⎤ ⎢ ⎥ M N s ⎢ ... ⎥ = 0 ⎢⎣ q& s ⎥⎦
]
(5)
where N p and N s are the primary and secondary network matrices, respectively, and q& p and q& s are the corresponding primary and secondary magnitude vectors, respectively. So, the secondary magnitude vector is given by: q& s = − N s−1 N p q& p
(6)
2.3. The Assur virtual kinematics chain concept The Assur virtual kinematic chain concept, virtual chain for short, is essentially a tool to obtain information about the movement of a kinematic chain or to impose movements on a kinematic chain (Campos et al. (2005)). This concept was first introduced by Campos (2004) which defines the virtual chain as a kinematic chain composed of links (virtual links) and joints (virtual joints) satisfying the following three properties: a) the virtual chain is open; b) it has joints whose normalized screws are linearly independent; and c) it does not change the mobility of the real kinematic chain, in other words, it is an Assur group (Baranov, 1985). The two useful virtual chains to describe movements in bi-dimensional space are the PPR and the RPR chain. 2.4.1 The orthogonal PPR Assur virtual chain The PPR virtual chain is composed by two virtual links (C1, C2) connected by two prismatic joints, whose movements are in the x and y orthogonal directions, and a rotational joint, whose the movement is in the z directions, see Fig.2. The prismatic joints are called px and py, and the rotative joint is called rz. The first prismatic joint (px) and the rotative joint (rz) are attached to the chain to be analyzed (real chain). Joint px connects the link R1 with virtual link C1, joint py connects virtual link C1 with virtual link C2, and joint rz connect the virtual link C3 with real link R2 (see Fig.2). Let the twist $ px represent the movement of link C1 in relation to link R1, twist $ py represent the movement of link C2 in relation to link C1, twist $ rz represent the movement of link R2 in relation to link C2. Therefore, the movement of link R2 in relation to real link R1 may be expressed by $ px + $ py + $ rz .
Figure 2. PPR Assur virtual chain. Consider the C-reference system (C-system) attached to the virtual link C2 at the rz joint. Therefore, there is no rotation between the C-system and the B-system (attached to the inertial base), and the rz joints are aligned with z axes. So, the normalized screws corresponding to the virtual joints represented in the C-system are
C
⎡1⎤ ⎡0 ⎤ ⎡0 ⎤ ˆ$ = ⎢0⎥ ; C $ˆ = ⎢1⎥ ; C $ˆ = ⎢0⎥ rz px py ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⎣0⎥⎦ ⎢⎣0⎥⎦ ⎢⎣1⎥⎦
(7)
It may be observed that the orthogonal PPR Assur virtual chain represents the movements in a planar Cartesian system. 2.4.2 The orthogonal RPR Assur virtual chain The RPR virtual chain is composed by one virtual links (C1,C2) connected by two rotatives joints, whose movements are in the z directions, and a prismatic joint, whose the movement is in the radial direction, see Fig.3. The prismatic joints is called pr, and the rotatives joints are called rz1 and rz2. The two rotatives joints (rz1 and rz2) and the prismatic joint (pr) are attached to the chain to be analyzed (real chain). Joint rz1 connect the link R1 with virtual link C1, joint pr connects virtual link C1 with virtual link C2, and joint rz2 connect the virtual link C3 with real link R2 (see Fig.3). Let the twist $ rz1 represent the movement of link C1 in relation to link R1, twist $ pr represent the movement of link C2 in relation to link C1, twist $ rz represent the movement of link R2 in relation to link C2. Therefore, the movement of link R2 in relation to real link R1 may be expressed by $ rz1 + $ pr + $rz 2 .
Figure 3. RPR Assur virtual chain. Consider the C-reference system (C-system) attached to the virtual link C2 at the rz2 joint. Therefore, there is a rotation between the C-system and the B-system (attached to the inertial base). So, the normalized screws corresponding to the virtual joints represented in the C-system are
C
⎡1 ⎤ ⎡0 ⎤ ⎡1 ⎤ ˆ$ = ⎢ 0 ⎥ ; C $ˆ = ⎢1⎥ ; C $ˆ = ⎢0⎥ pr rz 2 rz1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢q pr ⎥ ⎢ ⎢⎣0⎥⎦ ⎥ 0 ⎣ ⎦ ⎣ ⎦
(8)
where the qpr is the pr joint displacement . It may be observed that the orthogonal PPR Assur virtual chain represents the movements in a planar Cartesian system. Other Assur virtual chains can be found in Campos (2004) and Campos et al. (2005). 3. PROPOSED METHOD The proposed method consists in introduce an Assur virtual chain between the collision point of the manipulator and the obstacle, which will be used to monitor the obstacle approximation and to impose movement in order to avoid collision. To this end we define two different stages: monitoring and avoidance. In monitoring stage we intend to identify when the position of a joint (or joints) of the virtual chain is near to a limit given by the obstacle position. This joint position is obtained by the integration of differential kinematic where this virtual joint is obtained as a secondary joint. When the position of the virtual joint indicates the proximity of the obstacle the stage is changed to the avoidance stage. In the avoidance stage the joint used to monitor the collision possibility (as a secondary joint) is taken as a primary joint in order to impose a movement that avoids the collision. This changes the closed kinematic chain. To present the method, a 2D avoiding collision application is considered in the sequence. A 2D avoiding collision application Let a planar robot P3R, that represent a planar version of the robot of ROBOTURB project, composed by 5 links and four joins. The links are numerated from 1, at the base, to 5 on end-effector; and the joins are identified from A to D, see Fig.4. Let an obstacle defined by a point in the workspace as in Fig.4. The A joint is a prismatic joint and it movement is represented by the screw $A. The joins B, C and D are rotational joints and its movements are represented by the screws $B, $C, e $D, respectively.
Figure 4. Planar ROBOTURB Inverse differential kinematics without obstacle avoidance
Figure 5. Closed kinematic chain
To solve the inverse differential kinematics we introduce a PPR Assur virtual chain to close the kinematic chain (Fig.5) as in Campos et al. (2005), and apply the Davies method to determine the kinematic constraint equation (Eq. (4)),
[
N p = − $ px N s = [$ A
$ py $B
$C
$ rz
]
; q& p = [q& px q& py
$ D ] ; q& s = [q& A
q& rz
q& B
]
q& C
T
(9) T q& D ]
(10)
To solve the differential kinematics the secondary matrix has to be inverted (see Eq. (6)). In Eq. (10) it has more columns than rows, which characterizes a redundant manipulator. So, a method to solve redundant manipulators has to be used. In this example we use the method of kinematic constrains introduced in Santos et al. (2006). According the method of kinematic constraints we impose the movement of joint A and, consequently, take it as a primary joint instead secondary as before. So, we obtain a new kinematic constraint equation given by:
[
N p = − $ px
N s = [$ B
$ py
$C
$ rz
$A
]
; q& p = [q& px q& py
$ D ] ; q& s = [q& B
q& C
q& rz
q& A
]
T
T q& D ]
(11) (12)
Now, even the matrix N s is invertible, the secondary velocity vector q& s could be calculated using Eq. (6), Eq. (11) and Eq. (12) and giving the primary trajectory to be executed. Inverse differential kinematics with obstacle avoidance To describe the method proposed in this paper it is considered that we intend to avoid the collision of the point corresponding to the link 3 with the obstacle depicted in the Fig. 4. To this end, the first step is to introduce an Assur virtual chain connecting the manipulator’s collision point (in link 3) and the obstacle. As could be observed in Fig.6, the obstacle is considered as an extension of link 1. We chose a RPR virtual chain, described in section 2, because it allows to monitor or to impose the limit to the radial distance directly. As we wish to avoid the obstacle collision during a simultaneous trajectory execution, we should consider a virtual chain attached to the end effector, introduced to this end as before. The resulting kinematic chain is shown in Fig.6.
Figure 6. Planar robot PRRR connected with Assur virtual chains PPR and RPR (to collision avoidance) The kinematic chain shown in Fig.6 has two closed circuits. The robot chain and the PPR virtual chain used to define the end-effector desired movement compose the circuit 1. The partial robot chain until link 3 and the RPR virtual chain form the circuit 2, taken into account that the obstacle point is an extension of link 1. Applying the Davies method we determine the kinematic constraint equation (Eq. (4)) where:
⎡$ˆ $ˆ B N =⎢ A ⎢⎣$ˆ A $ˆ B
[
q& = q& A
q& B
$ˆC
$ˆ D
0
0
q&C
q& D
0 − $ˆ
rz1
0 − $ˆ
q& rz1 q& pr
pr
q& rz 2
0 − $ˆ
− $ˆrz
− $ˆ px
0
0
rz 2
q& rz
q& px
q& py
− $ˆ py ⎤ ⎥ 0 ⎥⎦ 6 ×10
(13)
]
T
(14)
Now, in monitoring stage, we choose joints px, py, rz and A as the primary ones as before. So, the resulting secondary joint velocity vector is calculated using Eq. (6), where:
[N ]
⎡$ˆ =⎢ A ⎢⎣$ˆ A
− $ˆ px
− $ˆ py
0
0
[N ]
⎡$ˆ = ⎢ˆ B ⎢⎣$ B
$ˆC
$ˆ D
0
0
p m
p m
− $ˆ rz ⎤ ⎥ ; q& p 0 ⎥⎦
[ ] = [q&
0 0 $ˆ rz1 $ˆ pr
q& px
q& py
q& rz
0 ⎤ & & ˆ$ ⎥⎥ ; [q s ]m = q B rz 2 ⎦
q& C
q& D
m
A
[
]
T
(15)
q& rz1
q& pr
q& rz 2
]
T
(16)
Note that the secondary matrix is square and even it is invertible the secondary vector [q& s ]m could be calculated using Eq. (6). So, the prismatic virtual joint velocity could be obtained. By integrating q& pr we obtain the prismatic joint position q pr , which is compared with a given limit, until it reaches this limit. Then the avoidance stage begins. In the avoidance stage, we choose the radial virtual joint pr as a primary one, in order to use it to impose the collision avoiding movement. At the same time, to maintain the resulting secondary matrix square, the robot joint A velocity is taken as a secondary one. So, joints px, py, rz and pr are the primary ones in the avoidance stage. Consequently, the resulting secondary joint vector is calculated using Eq. (6), in which:
[N ]
⎡$ˆ = ⎢ pr ⎢⎣ 0
− $ˆ px 0
[N ]
⎡$ˆ =⎢ B ⎢⎣$ˆ B
$ˆC 0
p a
p a
− $ˆ py 0
$ˆ D 0
− $ˆrz ⎤ ⎥ ; q& p 0 ⎥⎦
[ ] = [q& a
pr
q& px
0 $ˆ A 0 ⎤ ⎥ ; [q& s ]a = [q& B ˆ$ $ˆ A $ˆ rz 2 ⎥⎦ rz1
q& py
q& C
q& rz
q& D
]
T
q& rz1
(17)
q& A
T q& rz 2 ]
(18)
The secondary matrix [N s ]a is square and even it is invertible the secondary vector [q& s ]a could be calculated using Eq. (6). So, the whole robot joints movement, i.e., the joint velocity magnitudes q& A , q& B , q& C and q& D , is determined satisfying the trajectory to be executed and the collision avoiding requirements. 4. NUMERICAL SIMULATION To ilustrate the method we present a numerical simulation of the robot P3R working in a constraint environment defined by a punctual obstacle. 4.1 Simulation description Consider the planar robot P3R in Fig. 7 whose link lengths are link3 = link4 = 2.0m and link5 = 1.0m. The length of link2 is the displacement of joint A. The axis of prismatic joint A is rotated in the plane xy by an angle equal to θ / 4 rad about z axis. The initial posture is q0 = [1.36m 0.92rad − 1.33rad − 0.64rad ] , corresponding to the end-effector location: p = [4.5 0.9]T m, φ = [− π / 3] rad (see Fig.7).
Figure 7. Planar robot PRRR – Numerical simulation The motion trajectory is defined by a circular path with 0.5m radius and the center point on (4.5,0.9)m (Fig.7) specified by: ⎡ ⎛ πt ⎞ ⎤ ⎢ 4 + 0.5 cos⎜ 2 ⎟ ⎥ ; φ (t ) = ⎡− π − 0.4 sin ⎛ πt ⎞⎤ ⎜ ⎟⎥ ⎝ ⎠⎥ d ⎢ 3 pd (t ) = ⎢ ⎝ 2 ⎠⎦ ⎣ ⎢0.9 + 0.6 sin ⎛⎜ πt ⎞⎟⎥ ⎢⎣ ⎥ ⎝ 2 ⎠⎦
(19)
The collision region is defined by a circle path with 0.8m radius and the center point on (2.1,3.5)m; and the integration is performed in discrete time by resorting to numerical techniques, as depicted in (Guenther, 2007). 4.2 Graphical results Considering the aspects presented on section 4.1, Fig. 8 and Fig. 9 shows the results of inverse kinematic solution of robot P3R, based on the method proposed. On the monitoring stage (Fig. 8 and Fig. 9), the solution is obtained using the Eq. (15) and Eq. (16), and on the avoidance stage, the solution is obtained using the Eq. (17) and Eq. (18).
Figure 8. Inverse kinematic solution of P3R robot – Real joints.
Figure 9. Inverse kinematic solution of P3R robot – Joint Pr 5. CONCLUSION This paper presents a new methodology to solve the avoiding collision problem in robot manipulators. The proposed methodology is based on the method of kinematic constraint, of the screw theory, of Kirchhoff-Davies circulation law, and on the use of Assur virtual chains. The collision avoidance methodology was applied on a PRRR redundant manipulator, that’s the planar version of a 7 DOF manipulator, called ROBOTURB. The performed simulation clearly demonstrates the usefulness of the proposed methodology. 6. ACKNOWLEDGEMENTS This work was supported by by ‘Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPQ)’ Brazil. 7. BIBLIOGRAPHY Bottema, O. O., Roth, B., 1979, “Theoretical kinematics”, New York: North-Holand Pub. Co., pp. 533-544. ISBN 0444851240. Campos, A,“Cinemática Diferencial de Manipuladores Empregando Cadeias Virtuais”. 157 p. Tese(Doutorado em Engenharia Mecânica) — Universidade Federal de Santa Catarina, Florianópolis, Março 2004. Campos, A., Guenther, R. and Martins, D., 2005, “Differential kinematics of serial manipulators using virtual chains”, J. Brazilian Soc. Mechanical Sciences & Engineering, Vol. 27, No. 4, pp. 345-356. Caselli S., Reggiani, M. and Sbravati, R., ”Parallel path planning with multiple evasion strategies”. In: .Washington, DC: Proceedings of the IEEE Internacional Conference on Robotics and Automation, 2002. p. 260–266. Ceccarelli, M., 2000, “Screw axis defined by Giulio Mozzi in 1763 and early studies on helicoidal motion”, Mechanism and Machine Theory, Vol. 35, No. 6, pp. 761-770. Davidson, J. K., Hunt, K. H., 2004, “Robots and screw theory: applications of kinematics and statics to robotic”, New York: Oxford University Press Inc. ISBN 0-19-856245-4. Davies, T. H. ,“Kirchhoff’s circulation law applied to multi-loop kinematic chain”. Mechanism and Machine Theory, v. 16, p. 171–173, 1981. Guenther, R., Simas, H. and Pieri, E. R. ,“Concepção cinemática de um manipulador para volumes de trabalho restrito”, In: Congresso Nacional de Engenharia Mecânica, 2000, Natal/RN. 2000. Guenther, R., Cruz, D. F. M., Martins, D., Simas, H., 2007, “An integration method for differential inverse kinematic of closed chain robots”, Muller, A. ,“Collision avoiding continuation method for the inverse kinematics of redundant manipulators”. In: .New Orleans, LA: Proceedings of the IEEE Internacional Conference on Robotics & Automation, 2004. v. 2, p.1593– 1598. Piaggio, M. and Sgorbissa, A. ,“Ai-cart: an algorithm to incrementally calculate artificial potencial fields in real-time”. In: . Monterey, CA: Proceedings of the IEEE International Conference on Robotics and Automation, 1999.p. 238– 243.
Santos, Carlos H. F. dos, Guenther, Raul, Martins, Daniel et al, 2006 “Virtual kinematic chains to solve the underwater vehicle-manipulator systems redundancy”, J. Brazilian Soc. Mechanical Sciences & Engineering, Vol. 28, No. 3, pp. 354-361. Soucy, M., “Manipulator Path Planning with Multi-Resolution Potencial Fields and Fuzzy Logic Control”. Dissertação (Mestrado), University of Ottawa, March 2005. Soucy M. and Payeur, P. ,“Flexible fuzzy logic control for collision-free manipulator operation”. In: . Niagara Falls, Canada: Mechatronics and Automation, 2005 IEEE International Conference, 2005. v. 2, p. 723–728. Tian, L. and Mao, Z. ,“Fuzzy neuro controller for a two-link rigid-flexible manipulator system”. In: Singapore: Proceedings of the International Conference on Neural Information Processing, 2002. p. 1867–1871. Tsai, L., “Robot analysis: the mechanics of serial and parallel manipulators”, John-Wiley and Sons, New York, 1999. 8. RESPONSIBILITY NOTICE The authors are the only responsible for the printed material included in this paper.
