Bertrand Tondu
email: bertrand.tondu@insa-toulouse.fr

Serge Ippolito

Jérémie Guiochet

Alain Daidié

A Daidie

Keywords: anthropomorphic robot-arm, artificial muscle, McKibben muscle

A Seven-degrees-of-freedom Robot-arm Driven by Pneumatic Artificial Muscles for Humanoid Robots

Introduction

The recent development of robotic applications in nonindustrial fields, such as the current research in so-called humanoid robotics, needs to augment the friendly character of robot arms so that they interact with human subjects in total safety. A promising secure way of developing robot-arms adapted to proximity with humans consists of using artificial muscle actuators for which compliance is identical to that of human joints. If chemo-mechanical artificial muscles are still too far removed from practical applications [START_REF] Hebert | Dynamics and thermodynamics of artificial muscles based on nematic gels[END_REF][START_REF] Martin | Electrostriction in field-structured composites: basis for a fast artificial muscle[END_REF][START_REF] Bar-Cohen | Electroactive polymers as artificial muscles: a review[END_REF][START_REF] Kaneko | Polymer gels as soft and wet chemomechanical systems-an approach to artificial muscles[END_REF], pneumatic artificial muscles, notably braided ones, for which the so-called McKibben artificial muscle is the most interesting representative, have proved efficient [START_REF] Schulte | The characteristics of the McKibben artificial muscle. The Application of External Power in Prosthetics and Orthotics[END_REF][START_REF] Matsushita | Synthesis of rubber artificial muscle[END_REF][START_REF]Rubber muscles take robotics one step further[END_REF][START_REF] Inoue | Rubbertuators and applications for robots[END_REF][START_REF] Caldwell | Control of pneumatic muscle actuators[END_REF]Tondu andLopez 1995, 2000;[START_REF] Chou | Measurement and modeling of McKibben pneumatic artificial muscles[END_REF][START_REF] Caldwell | Bio-mimetic principles in actuator design for a humanoid robot[END_REF][START_REF] Davis | Enhanced modeling and performance in braided pneumatic muscle actuators[END_REF]. Current industrial pneumatic artificial muscles are essentially derived from the McKibben model (Inoue 1968; see Fluidic Muscle MAS, technical document of Festo AG & Co. KG, Esslingen, Germany; see also The SHADOW Air Muscle, Shadow Robot Group, http://www.shadow.org.uk). With its high power-to-weight and power-to-volume ratios, closely imitating the functions of the natural skeletal muscle, the McKibben muscle is the most adapted artificial muscle for motorizing "soft" and human-size robot arms. Tire manufacturer Bridgestone has shown the possibility of designing and controlling robot-arms entirely actuated by McKibbentype artificial muscles called "rubbertuators" (i.e., rubber actuators) in which revolute joints are driven by two antagonistic muscles. Bridgestone "soft arms" are generally limited to two degrees of freedom (DoF; Noritsugu, Tanaka, and Yamanaha 1996), four DoF (Bridgestone Corporation 1987;[START_REF] Pack | A rubbertuator-based structure-climbing inspection robot[END_REF] or five DoF (Bridge-stone Corporation and Taicubo Engineering 1993). Other robots equipped with McKibben-type artificial muscles are structures limited to a few DoF, such as the ISAC Vanderbilt humanoid robot formed of two 3-DoF arms [START_REF] Kawamura | Intelligent robotic systems in service of the disabled[END_REF], the Salford bipedal system composed of two 3-DoF legs and a 2-DoF hip [START_REF] Caldwell | Investigation of bipedal robot locomotion using pneumatic muscle actuators[END_REF], or limited to a 1-DoF joint robot in the abstract [START_REF] Tuijhof | Design, actuation and control of an anthropomorphic robot arm[END_REF]. Our aim is to study the feasibility of designing a 7-DoF anthropomorphic robot-arm motorized by McKibben muscles, within this dual industrial and academic framework. Very recently, the Salford University Robotics team published results concerning the development of a 7-DoF exoskeleton aimed at the rehabilitation and training of an upper limb, motorized by pneumatic McKibben muscles [START_REF] Tsagarakis | Development and control of a "soft-actuated" exoskeleton for use in physiotherapy and training[END_REF]. This exoskeleton, which is lighter than 2 kg, has been successful, but resembles more an orthotic device than a robot-arm. On the other hand, it is strictly speaking a 6-DoF structure due to a joint coupling at the shoulder, whereas our 7-DoF anthropomorphic arm ensures true independence between shoulder abductionadduction and flexion-extension motions. In Section 2 of this paper we define the anthropomorphic architecture, in Section 3 we introduce the artificial muscle actuator used, in Section 4 we detail the mechanical design of the robot, and in Section 5 we demonstrate the testing of the teleoperated arm using a pair of joysticks.

Anthropomorphic 7-DoF Robot Kinematic Architecture

An industrial robot, the tool of which is to be positioned and oriented in space, is typically a six DoF articulated chain. Considering the robot's working area and dexterity, one of the most interesting architectures of the six DoF is the purely revolute architecture associating a 3R carrier known as "revolute type" [START_REF] Rosheim | Robot Evolution. The Development of Anthrobotics[END_REF] or "anthropomorphic type" with a 3R wrist structure essentially of the roll-pitch-roll type. Industrial RX robots derived from classical PUMA robots follow this model. Their wide joint ranges generate a working area that can be globally as wide as that covered by the two arms of a human operator. However, this universal robot has, as a kinematic fault, position singularities specific to the shoulder and elbow, as well as an orientation singularity specific to the wrist [START_REF] Hollerbach | Optimum kinematic design for a seven degrees of freedom manipulator[END_REF]. These singularities reduce the robot's ability when Cartesian trajectories have to be performed. The interest for researchers in robotics of redundant structures arises from the possibility of both solving the question of singularity and avoiding obstacles [START_REF] Kreutz-Delgado | Kinematic analysis of 7-DOF manipulators[END_REF]. Difficulties for controlling redundant robot-arms lead to preferring 7-DoF revolute structures, which are here called 7R robot-arms. One of the 7R's most interesting structures is that derived from the anthropomorphic 6R structure by adding an elbow roll. The JPL "dexterous arm" [START_REF] Seraji | Motion control of 7-DOF arms: the configuration control approach[END_REF] and the indus-trial Mitsubishi-PA 10 arm were designed according to this structure. The seven DoF of this structure can be described in roll-and-pitch terms by distinguishing three anthropomorphic joint sets performing the following motions: shoulder roll and pitch, elbow roll and pitch, and wrist roll and pitch associated to the tool plate roll.

Furthermore, this kinematic structure can be interpreted as a simple model of the kinematic structure of the human arm, for which the shoulder roll-and-ball socket is modeled by the first three DoF, and the wrist roll-and-ball socket the last three DoF Within the framework of this paper, we propose to apply the terminology of the joint physiology to this anthropomorphic 7R structure rather than the roll-pitch-yaw mechanical terminology. Joint physiology specifies the motions of body links in terms of external-internal rotation, flexion-extension, and abduction-adduction. The use of this terminology implies the choice for a reference body position, called the zeroanatomical position, which is defined as follows. The person stands in an upright position, arms resting along the body, palms turned to the front. In this position, the spatial rotation motions of a mobile link are defined as follows: externalinternal rotation is the rotation of the link along its long axis; flexion-extension is the motion that decreases-increases the angle between jointed links; abduction-adduction moves the link away-towards the body sagittal plane (i.e., dividing the body between its right and its left). Figure 1 illustrates the correspondence which can be established between roll-pitchyaw notation and joint physiology motion notation.

This physiological approach to kinematic robot-arm structures is suggested by the development of humanoid robots in which arms have a first horizontal axis performing flexion-extension or abduction-adduction movements instead of the central vertical axis of industrial robot-arms. For example, the arms of the Honda humanoid robots P2, P3 or Asimo have a 7R structure for which the first axis performs a shoulder flexion-extension [START_REF] Hirai | The development of Honda humanoid robot[END_REF]; http://world.honda.com/ASIMO/). The choice for the first axis-whether shoulder flexion-extension or shoulder abduction-adduction-is free. However, in order to be in agreement with the joint physiology definition of the horizontal flexion-extension [START_REF] Kapandji | Physiology of the Joints[END_REF], it was decided to consider the shoulder abduction-adduction first. Consequently, if abduction is zero, shoulder flexion-extension is performed on a vertical plane, and on a horizontal plane if abduction is 90 • , as illustrated in Figure 2.

Finally, the following sequence was considered for the shoulder: abduction-adduction, flexion-extension, and external-internal rotation. For reasons of symmetry, first abduction-adduction and secondly flexion-extension were also considered for the wrist. Figure 3 shows this 7R structure, where it can be noted that the wrist corresponds in mechanical terms to a roll-pitch-yaw wrist type instead of the roll-pitch-roll wrist type of usual 7R anthropomorphic robotarms. Figure 3(a) specifies the association of the rotation axis with the bodily structure considered in the form of the three main segments: arm, forearm, and hand, in the so-called zeroanatomical position where the arm falls along the body with the palm open to the front. The positive sense of every joint axis was chosen to correspond to an "opening" motion of the upper limb, as specified in Table 1. Because the wrist is fully stretched in the zero-anatomical position, a zero-robot position is considered in which the wrist joints are in a medium configuration, as illustrated in Figure 3(b). The robot hand in use is a simple two-finger grip mimicking the thumb in opposition to the set of other fingers. Before dealing with the mechanical design of the anthropomorphic robot-arm, the main properties of the McKibben muscle actuators are described, which are to motorize this robot-arm.

Naturally Compliant Artificial Muscle Actuators

Fundamental Model of Artificial Muscle

Skeletal muscles are characterized by a contraction force depending on the length of the muscle. During physiological contraction, muscle tension decreases in the same way as muscle length. Consequently, the skeletal muscle behaves like a spring for which the stiffness depends on the nervous control of the muscle. Although the nervous system can vary the muscle force in a complex manner, either by recruitment of motor units or modulation of the firing rate of motor neurons [START_REF] Ghez | Muscles: effectors of the motor system[END_REF], the following fundamental linear model of the skeletal

F = uf max (1- ε ε max ), 0 ≤ ε ≤ ε max and 0 ≤ u ≤ u max ε = (l 0 -l)/ l 0 .
(1)

Here, u designates the value of the nervous control on a [0, u max ] range, ε is the contraction ratio of muscle length l with regard to the resting muscle length l 0 in a [0, ε max ] range, and f max is the maximum force per control unit. From eq. ( 1), the muscle stiffness equation is obtained:

F r = -u f max ε max δε. ( 2
)
This means that the muscle, deviating from its equilibrium position of a δε, returns to this position due to a restoring force F r , as illustrated in Figure 4. This model can be used as a reference model for any type of artificial muscle which is defined as a spring-like device for which the control variable is proportional to its stiffness. A competitive approach defines the notion of artificial muscle from the skeletal muscle in Hill's model, which favors a tension-velocity relationship (Klute, Czerniecki, andHannaford 1999, 2002). We preferred a tension-length model so as to highlight the "natural compliance" of the artificial muscle actuator resulting from the muscle's spring-like behavior. On the other hand, as opposed to a Hill-based model, using a spring-type model more easily distinguishes the stiffness component from the viscous-type component, which will be considered later. Any artificial muscle must be naturally damped but, especially in the case of fluid artificial muscles, this natural damping depends on the materials used to build the artificial muscle (inner tube, textile braid) as well as the servo-valve supplying the muscle. In particular, the McKibben muscle, as a system where input is the control pressure and output the muscle length, can appear as poorly damped [START_REF] Klute | Artificial muscles: actuators for biorobotic systems[END_REF]. This means that no damping would occur in response to a physical pressure step. However, in practice, muscle pressure is generated by using a specific pneumatic or hydraulic servo-valve contributing, in conjunction with a judicious choice of textile braid, to naturally damping dynamic muscle behavior (see Figure 5(d)).

McKibben Artificial Muscle

In the fundamental length-tension eq. ( 1) model, no hypothesis is made on the nature of control variable u. The notion of the artificial muscle is consequently very free concerning the physical principle that produces the contraction force. However, in previous studies we have shown how beneficial it is to use air-pressurized inner tubes surrounded by braided shell to design artificial muscles adapted to robot-arms. The socalled McKibben artificial muscle today appears as one of the most promising artificial muscles, as proved by the numerous papers devoted to it. Basically, when the McKibben muscle inner tube is under a pressure P , the double-helix braided shell surrounding it transforms the circumferential stress of the inner tube into an axial contraction force F , which can be expressed in the fundamental generator force model [START_REF] Tondu | Modeling and control of McKibben artificial muscle robot actuators[END_REF]:

F (ε, P ) = (π r 2 0 )P [a(1 -kε) 2 -b], 0 ≤ ε ≤ ε max ε = (l 0 -l)/ l 0 , a = 3/ tan 2 (α 0 ), b = 1/ sin 2 (α 0 ). (3)
Instead of contraction ratio ε, an equivalent model considers the current braid angle (i.e., the angle between the muscle axis and one strand of the textile weave) as the muscle characteristic contraction variable [START_REF] Schulte | The characteristics of the McKibben artificial muscle. The Application of External Power in Prosthetics and Orthotics[END_REF][START_REF] Chou | Measurement and modeling of McKibben pneumatic artificial muscles[END_REF][START_REF] Davis | Enhanced modeling and performance in braided pneumatic muscle actuators[END_REF]. It seems to us preferable to favor contraction ratio ε (or the current muscle length) as an output variable in order to more easily derive the antagonistic muscle actuator model. The proposed model characterizes McKibben muscle of current length l by its inner tube, initially cylindrical, of initial length l 0 and of initial radius r 0 , and by its textile weave of initial braid angle α 0 , which is derived from the data of our cotton braid supplier. 1 Note that the use of cotton gives a natural softness to the braid, which helps to completely cover the inner tube with the braid at rest state. In these conditions, McKibben artificial muscle is clearly geometrically characterized in its non-stretching rest state by the three parameters l 0 , r 0 , and α 0 . Conversely, the use of plastic-type fibers, as considered by other researchers, creates a braid stiffness which "opens" the braid and makes it stretchable from its rest state to the state in which the braid completely covers the inner tube. There results a real difficulty in equating muscle initial state l 0 , r 0 with initial braid angle α 0 as considered here.

However, the purely cylindrical model with three parameters is too simple to express the artificial muscle behavior as a force generator in a satisfactory manner. This is why we have introduced a correction parameter k. Factor k is an empirical factor, experimentally derived, expressing both a bound effect (muscle does not maintain its cylindrical shape at end proximity when it contracts) and the non-integral transmission of pressure force into the muscle weave. Modeling muscle envelope evolution by finite elements, as envisaged by the Washington Biorobotics Laboratory (http://brl.ee.washington.edu), is a promising way to theoretically model the bound effects, but might not be sufficient to fully explain the high dependence between control pressure and real maximum contraction ratio. Faced with this complexity, we preferred to preserve a simple shape of the muscle force generator model by introducing the empirical parameter k, which can be estimated as a constant or as a function dependent on pressure P .

Compared to the general artificial muscle eq. ( 1) model, the McKibben muscle tension-length model appears not to be as exactly linear as eq. ( 3) expresses it. It does however attain the desired artificial muscle nature by generating stiffness globally proportional to the control pressure. We have developed our own artificial McKibben muscles and tested them on the experimental setup, which are used to record muscle parameters during isometric and isotonic contraction as defined by physiology [START_REF] Ghez | Muscles: effectors of the motor system[END_REF]. Figure 5(a) shows an isometric contraction for a muscle of 120 mm initial length, 7 mm initial radius, and a 23 • initial braid angle, aimed at motorizing the wrist of our 7-DoF prototype. Figure 5(b) shows isometric contraction for a muscle of 230 mm initial length, 12 mm initial radius, and a 17 • initial braid angle, aimed at motorizing the robot-arm shoulder. The resulting experimental static tension-length at constant pressure emphasizes the ability of McKibben muscle to develop high maximum force in relation to its own weight and volume. Under 5 bars, the considered muscle can develop a maximum force greater than 5000 N for 1. The application of textile technology classical formulae indeed gives tan(α 0 ) = π(D int -2d)/p, where D int is the internal diameter of the braid, d is the diameter of the thread, and p is the braid pitch that can be determined from p = n s /2n m (in cm), where n s is the number of spindles and n m is the number of meshes by centimetre.

a weight of about 100 g and a maximum volume of about 500 cm 3 . Note that the McKibben muscle maximum force can be dimensioned by imposing the value either of the initial braid angle or of the initial muscle radius. However, as is clearly shown by eq. ( 3), maximum muscle force is very sensitive to a small variation in the initial braid angle. Consequently, we preferred this parameter to define the set of muscles motorizing the robot, thus obtaining a wide range of maximum contractile force. It is important however to emphasize that the smaller the initial braid angle, the more the muscle radius increases during contraction, which thus requires an adequate choice of rubber able to tolerate the corresponding elongation.

However, the power of the McKibben muscle is limited by its pneumatic nature and particularly by the dynamic behavior of the servo-valve feeding it with pressurized air. The contraction of the McKibben muscle does indeed go with air consumption all the more higher as the initial braid angle is low. Considering a purely cylindrical model of the McKibben muscle (i.e., r = (r 0 / sin(α 0 )) √ 1cos2 (α 0 )(1ε) 2 ; Tondu 1995) the following approached expression of the current muscle volume V (ε) can be derived

V (ε) = V 0 (1 -ε)(1 -cos 2 (α 0 )(1 -ε) 2 )/ sin 2 (α 0 )
where V 0 is the muscle initial volume (π r 2 0 l 0 ). The result is that the ratio between the current muscle volume and its initial volume has a maximum value corresponding to the theoretical maximum contraction ratio 2 ε max = 1 -(1/ √ 3 cos(α 0 )) and equal to (2/3 √ 3 cos(α 0 ) sin 2 (α 0 )). The classic formula of perfect gases P V = (m/M)RT (where P is pressure, V is volume, m is mass, M is molar mass, R is gas constant, and T is temperature) can be applied to determine air consumption m between the zero-contraction state and the maximum contraction state for a P pressure variation. We obtain m = (2/3 √ 3 cos α 0 sin 2 α 0 )(V 0 M/RT ) P . In the case of the biggest shoulder muscles of our 7R robot (joint 2: l 0 = 230 mm, r 0 = 12 mm, α 0 = 17 • ) the muscle volume initially equal to about 104 cm 3 can theoretically increase until about 470 cm 3 , which for a P =5 bar pressure variation corresponds to an air consumption of 2.85 g, a little less however in practice. In the case of the servo-valves feeding the muscles of our robot-prototype-Samson I/P converters (Samson Corporation, Frankfurt, Electropneumatic Converter I/P 5288), generally considered to be rapid in their 0-5 bar output pressure range-a relatively slow pressure rise results, as shown in the example in Figure 5(c), corresponding to an isotonic contraction of the above-mentioned shoulder muscle picking up a 15 kg mass in response to a 5 bar numerical step. Corresponding muscle length evolution is given in Figure 5(d). During contraction, maximum muscle power can be estimated at 75 W, which appears as a relatively weak value given the maximum muscle force (in comparison, the skeletal 1987). The increasing of the McKibben muscle power necessitates, in consequence, the increasing of the bandwidth of the servo-valves feeding it. 3 Interesting proposals have been made with this end in view [START_REF] Davis | Enhanced modeling and performance in braided pneumatic muscle actuators[END_REF]. It is important however that the sizes of the valves, just as the sizes of the air pipes, do not become out of proportion to the muscle size.

The low bandwidth of the McKibben muscle corresponds to moderate response times. The McKibben muscle response 3. We have analyzed elsewhere the bandwidth of our I/P converters [START_REF] Boitier | Design and control of a two DoF SCARAtype robot actuated with McKibben pneumatic artificial muscles[END_REF]. Its value appears to be independent of the mean working pressure if it is greater than 2 bars but, naturally, depends on output volume; it can be estimated at 5 Hz for a 75 cm 3 volume, at 1 Hz for a 300 cm 3 volume, and at 0.5 Hz for an 800 cm 3 volume. Consequently, the bandwidth of our robot muscles varies as a function of the muscle size between about 1 and 4 Hz. time depends on multiple factors (e.g., initial muscle volume, controlled pressure variation, and load), but is generally between 0.2 and 1 s (Figure 5(d) gives a typical result for the strongest muscles of our robot which are also the slowest) when the skeletal muscle contracts between 10 and 300 ms [START_REF] Boff | Engineering Data Compendium: Human Perception and Performance, Harry G[END_REF]. This relative slowness of the artificial muscle contraction could appear as a real handicap for the development of robot-arms actuated with McKibben artificial muscles. However, it is important to note that artificial muscle robots aim at performing tasks at "human" speed and accuracy, and today's dynamic McKibben muscle performances are relatively close to those of their human skeletal equivalents. However, researchers have become interested in the abilities of the human arm to perform high-frequency tasks such as drum rolls at a 30 Hz frequency despite the relative slow response of the neuromuscular system [START_REF] Hajian | Drum roll: increasing bandwidth through passive impedance modulation[END_REF]. By testing their hypothesis on McKibben muscles, they have shown how the modulation of muscle passive impedance can explain this paradox. Such an approach could be applied to the control of McKibben artificial muscle robot actuators.

Furthermore, all researchers who have worked on braided muscles have emphasized the difficulty of obtaining an accurate model of muscle contraction (see, for example, [START_REF] Klute | Artificial muscles: actuators for biorobotic systems[END_REF]. This difficulty is mainly due to the structure itself of the artificial muscle consisting of soft materials whose interactions are complex, and also to the necessity of separating the specific muscle dynamics from the servo-valve dynamics. In a general way, McKibben muscle dynamic force can be written as

F dyn = F -F elastic -F f riction , ( 4
)
where F is the static force generated by the muscle as expressed in eq. ( 3), F elastic is an elastic force component associated with the inner tube dynamic behavior, and F f riction is a dry and kinetic complex friction force component associated with the braid dynamic behavior. Friction forces seem to be particularly important in explaining the specific actuator hysteresis [START_REF] Tondu | Modeling and control of McKibben artificial muscle robot actuators[END_REF]. In comparison with classical actuators, the dynamic complexity of the McKibben muscle appears to be the price to be paid for the natural softness conferred to the robot by the artificial muscle. It implies an increased but non-insuperable control difficulty of the artificial muscle antagonistic actuator.

Antagonistic McKibben Artificial Muscle Actuator

Based on the model of the natural biceps-triceps system, two artificial muscles can be set into antagonism for defining a revolute actuator according to the basic scheme of Figure 6(a). The working principle of the actuator is as follows. The two identical muscles are initially contracted of a ratio ε 0 = ε max /2 under a pressure of P 0 = P max /2, where P max is generally chosen equal to 5 or 6 bar and ε max is generally equal to a value between 0.1 and 0.15, due to a mean practical maximum contraction range of between 20% and 30% (see Figures 5(a) and (b)). From this initial situation, each muscle can contract or lengthen by ε 0 2. If the actuator drive is performed by a chain wheel of radius R, the motor torque generated by the actuator is given by

T = R[F 1 (ε 1 , P 1 ) -F 2 (ε 2 , P 2 )]. (5)
In this expression, F 1 designates the force produced by muscle 1 (defined as that producing positive actuator rotation motion) inflated at pressure P 1 and contracted of a ratio ε 1 = ε 0 -(Rθ/ l 0 ), and F 2 designates the force produced by muscle 2, inflated at pressure P 2 and contracted by a ratio ε 2 = ε 0 -(Rθ/ l 0 ). According to whether pressure P 1 > P 0 or P 2 > P 0 , F 1 or F 2 is an agonistic-type force, the force produced by the other muscle is an antagonistic-type force. From the force generator model of eq. ( 3), with a k constant typically between 1.25 and 1.35, the following expression of the static motor torque can be derived:

     T = K 1 (P 1 -P 2 ) -K 2 (P 1 + P 2 )θ , with K 1 = (π r 2 0 )R[a(1 -kε 0 ) 2 -b] K 2 = (π r 2 0 )R2a(1 -kε 0 )kR/ l 0 . ( 6
)
This model is particularly interesting as it shows how the antagonistic McKibben artificial muscle actuator performs the principle of the equilibrium point found in physiology [START_REF] Bizzi | Does the nervous system use equilibrium-point control to guide single and multiple joint movements?[END_REF]).According to this principle, the static motor torque can be divided into two effects-an agonistic effect T ago and an antagonistic effect T antago -as follows

T = T ago -T antago , with T ago = P ago (K 1 -K 2 |θ |) T antago = P antago (K 1 + K 2 |θ |) , (7)
where (P ago , P antago ) is equal to (P 1 , P 2 ) if θ ≥ 0 and to (P 2 , P 1 ) if θ ≤ 0. Equilibrium point θ equ is defined by equality T ago = T antago , which leads to the following expression of θ equ :

θ equ = K 1 (P 1 -P 2 ) K 2 (P 1 + P 2 ) , (8)
as Figure 6(b) illustrates. When the actuator deviates from the equilibrium point of value δθ, it returns to it by means of a restoring torque

T r = -K 2 (P 1 + P 2 )δθ, (9)
which expresses the actuator stiffness theoretically estimated by the term K 2 (P 1 +P 2 ). Generally speaking, the antagonistic McKibben muscle actuator is a system with a double input (P 1 , P 2 ) and a single output θ , for which the control would need a multivariable approach. However, it is possible to control the actuator as a single input-single output by defining agonistic and antagonistic muscle controls symmetrically as follows: P 1 = P 0 + P , P 2 = P 0 -P . Consequently, the static motor torque expression becomes a linear function of the input variable P and the output variable θ :

   T = k 1 P -k 2 θ,
where

k 1 = 2K 1 , k 2 = 2K 2 P 0 . (10)
It is important, however, to note that the muscle volume does not change linearly with its contraction ratio. As a consequence, the air consumption of the two antagonistic muscles is not symmetrical. 4 The result is a slight delay in the buildup of the desired pressure in the two antagonistic muscles, all the higher as the two muscles have different contraction ratios, i.e., actuator angle θ is close to its extreme values.

In the equilibrium position, the actuator behaves like a spiral spring whose stiffness is now independent of θ and equal to k 2 . This stiffness is particularly adapted to human contact: for example, the elbow flexion-extension joint 4 of our 7R robot is motorized by muscles of parameters l 0 = 195 mm, r 0 = 8.5 mm, α 0 = 19 • , k = 1.25, with a chain wheel of radius R = 14.5 mm and P 0 = 2.5 bar, ε 0 = 0.1, leading to an estimated stiffness of about 15 N.m.rd -1 near to values given by joint physiology [START_REF] Boff | Engineering Data Compendium: Human Perception and Performance, Harry G[END_REF]. The term "natural compliance" has been given to this joint compliance, which has the same nature as that of the joints of our skeleton. However, this natural compliance would be of no practical use if the actuator were not naturally damped. As previously emphasized, this natural damping of the actuator is due to particularly complex elastic and friction phenomena specific to the muscle components as it is to the servo-valve dynamic behavior. Despite this complexity (and apart from hysteresis behavior) the dynamic expression of the motor torque can be satisfactorily approached by a linear expression as follows

T dyn = a 1 P -a 2 θ -a 3 θ, (11)
leading to a second-order linear model of the actuator that drives a constant inertial load. Defined in this way, the McKibben artificial muscle actuator is an open-loop stable system. In consequence, the system is easily identifiable and the corresponding identified model can be used as a forward action of a linear or non-linear corrector [START_REF] Tondu | Modeling and control of McKibben artificial muscle robot actuators[END_REF].

On the one hand, the antagonistic McKibben artificial muscle actuator offers static performances in maximum force-muscle weight and maximum force-muscle volume ratios well adapted to the motorization of human-size robot-arms. It also offers dynamic performances in response time, damping and compliance relatively well adapted to generating humanlike motions, and, moreover, improvable. However, the maximal contraction that varies between 15% and 40%, depending on muscle parameters and control pressure, is a major disadvantage of the McKibben muscle. In comparison, the skeletal muscle contracts only 30% at rest, but in vivo the muscles of vertebrates are stretched by about 10-20% of their length at rest [START_REF] Boff | Engineering Data Compendium: Human Perception and Performance, Harry G[END_REF]. In vivo skeletal muscle can then contract by 40-50%. This structural limitation specific to McKibben muscles imposes, for a desired joint range, the obligation to look for an adequate compromise between muscle length and pulley radius with the possible help of a speed increaser, as will be seen in Section 5. Unfortunately, natural muscle compliance implies that motor torque is highly dependent on the joint variable. An artificial muscle robot is consequently highly dependent on gravity, unlike classical industrial robots, which are free from any gravity effect due to their speed reducers. The design of our anthropomorphic artificial muscle robot-arm aims at taking advantages of the McKibben muscle actuator properties.

Design of the 7R Anthropomorphic Robot-arm

The classic, and essentially electric, revolute actuators of industrial robotics are characterized as to bulk by a compact structure globally symmetrical around the actuator rotation axis. In comparison, the antagonistic McKibben muscle actuator shows very different bulk mechanical characteristics, as follows.

1. First, muscle length is a specific bulk factor because the actuator joint range directly depends on it; this factor is all the more significant as the muscle maximum contraction ratio is limited. Like the human muscularskeletal system, it is possible however to integrate the artificial muscles along the links of the robotized structure.

2. Secondly, the rotation axis is perpendicular to the actuator two-muscle "plane" and located at an edge of the pair of muscles. In consequence, the artificial muscle actuator does not show this axial symmetry peculiar to most revolute actuators.

As a consequence, the design of a robot actuated with artificial muscles must be adapted to the nonclassical character of this revolute actuator, in order to find the right compromise between acceptable dimensions of the robot and desired dynamic performances. The goal of our study consisted of designing a robot-arm where the arm and forearm segments would have lengths similar to those given by human biometry, and the joint ranges are defined according to the 7-DoF model of the human arm presented in Section 2, while minimizing the total weight of the robot. This minimization combined with the joint natural compliance is indeed an essential element for the safe functioning of the robot in a human environment.

Dimensioning of the McKibben Artificial Muscle Actuator

If we consider, as illustrated in Figure 6(a), that each muscle is initially contracted by ε 0 the actuator joint range is given by

θ range = [-ε 0 l 0 /R, + ε 0 l 0 /R] (rd). ( 12
)
Furthermore, in the functioning conditions given in Section 3.3, it appears that the maximum torque in absolute value is generated for a variable control P to P 0 and for an angular variable θ equal in absolute value to |ε 0 l 0 /R|. Derived from eq. ( 10), the following expression of the maximum torque results

T max = k 1 P 0 + k 2 (ε 0 l 0 /R),
which can be expressed in the following form by using eq. ( 6) 5

T max = 2(πr 2 0 )R[a(1 -k 2 ε 2 0 ) -b]P 0 . ( 13
)
Comparing eqs. ( 12) and ( 13), it appears that, if joint range and maximum torque naturally depend both on muscle parameters and driving pulley radius, the joint range does not 5. This result corrects the incomplete analysis of maximum torque generated by the McKibben muscle actuator we have previously developed [START_REF] Tondu | Modeling and control of McKibben artificial muscle robot actuators[END_REF].

depend on α 0 and r 0 Consequently, these latter two parameters can be specifically used for the dimensioning of the maximum torque. It is also important to note that the joint range directly depends on initial muscle length l 0 , but a lengthening of the muscles unfortunately increases the bulkiness of corresponding mechanical parts supporting the muscles. If pulley radius R plays the classic role of increasing or decreasing joint range and maximum torque in an inverse ratio, initial muscle contraction ε 0 plays a specific role this study would like to emphasize. Its value has a weak influence on the maximum torque and a large influence on the joint range.

As previously stated, a typical range for ε 0 is 0.1-0.15 but it can be interesting to adjust this value accurately, either to increase the maximum joint range by choosing ε 0 near the estimated maximum value, or to use ε 0 as a factor for tuning the corresponding joint mechanical stop. This approach will be used in our robot prototype. Furthermore, the current actuator torque appears to depend on the angular variable generated by the actuator, in analogy with human gestures. For the robot actuated with artificial muscles, a load manipulation results fundamentally different from that of industrial robots, which is generally performed at constant torque. The consequence of this joint dependence of the actuator torque will be analyzed in Section 5 when the robot moves against gravity.

Mechanical Arm Design Based on the Use of CAD Modeling 6

Mechanical arm design aims to develop technological solutions adapted to the artificial muscle actuator, which combines arm compactness globally equivalent to that of the human arm, with joint ranges approaching those of the human arm, as indicated in Figure 3(b).

In comparison with Bridgestone soft-arms and derived robots, the relatively advanced anthropomorphic character of our robot-arm has led us to develop original solutions for the shoulder and the wrist design in order to achieve specific abduction-adduction and flexion-extension motions. The θ 1 and θ 2 shoulder joints impose the most difficult muscle dimension constraints because the corresponding joint actuators have to develop the highest motor torques for driving the arm against gravity. Making these two joints independent is also a major difficulty, as already emphasized by [START_REF] Tsagarakis | Development and control of a "soft-actuated" exoskeleton for use in physiotherapy and training[END_REF] whose McKibben muscle exoskeleton shows a shoulder flexion-extension and abduction-adduction dependence. The original solution we have considered to obtain this independence is based on the possibility of contracting the two ends of the McKibben muscle simultaneously. Figure 7 illustrates the joint 1 and joint 2 corresponding mechanical system. The joint 1 muscle pair is fixed to the robot frame. It 6. Details of the mechanical design developed on I-DEA CAD software are not given within the framework of this paper. A more specifically mechanical paper presenting the design of our anthropomorphic arm is being prepared. controls a pair of spur gears (ϕ and κ), driven by the actuator pulley λ, which amplifies the actuator angular motion in a 1.5 ratio. In opposition, joint 2 muscles can contract their two ends simultaneously. Joint 2 moves according to a differential principle performed by two bevel gears, functioning in opposition by means of a back pulley mechanism. When joint 2 muscles are controlled, they both act on pulleys µ and ν generating corresponding opposite angular displacements (-θ 2 and +θ 2 ). Consequently, joint 2 can move independently of joint 1, and when joint 2 does not move, the joint 2 driving system acts as an extensible tie that prevents pulleys µ and ν from rotating. As in the case of joint 1, the joint 2 bevel gear system amplifies the actuator angular motion in a 1.5 ratio. It is however important to note that if the considered mechanical system of joint 2 can generate a torque equivalent to that which would be obtained with a fixed end, the corresponding joint range is divided by 2 since it corresponds to twice the muscle contraction when one end is fixed. This joint range deficiency could be compensated by doubling the muscle length but would lead to unduly increasing the robot frame length. Consequently, joint 2 muscles are very slightly oversized compared to joint 1 muscles with a corresponding ε 0 tuned to its maximum value. As mentioned in Table 2, we can therefore hope to expect a 180 • joint 1 range and a 100 • joint 2 range with maximum torques greater than 60 N m.

Wrist flexion-extension and abduction-adduction are not so severely constrained in torque and joint range performances. Joint 6 is directly driven by its muscle pair as illustrated in Figure 8(a), and joint 7 is independently directly driven by means of a swing-bar of half length L 7 size which plays the role of the pulley radius R of other joint actuators, as illustrated in Figures 8(b) and(c). This mechanical solution gives the wrist a compactness in accordance with the desired human-like robot dimensions.

Table 2 synthesizes the set of muscles and transmission system characteristics. They were selected so as to generate anatomically desired joint ranges with associated maximum possible torque, while maintaining robot-arm dimensions similar to those of human arms. Note that the given value of "peak force" can be estimated from torque eq. ( 13) (with a typical k equal to 1.25) divided by the ratio m V . Figure 9 gives the resulting prototype. Its arm, in the anatomical sense, has a 351 mm length between axes, and its forearm a 307 mm length between axes. The lateral frame supporting the shoulder muscles is 350 mm wide. The hand is a simple two-finger grip actuated by a double-effect pneumatic cylinder. Figure 10 shows the current development of the robot with its muscle actuators. The total weight of the robot's mobile parts including muscles and chains is about 6 kg and the robot's total weight with its shoulder (fixed base) is about 10 kg. The weight of the pneumatic interface, which is essentially the weight of the 14 I/P converters supplying the muscles as the I/P converter supplying the robot tool-plate, is about 10 kg. These converters have been set in a lateral box attached to the robot frame. The size of used I/P converters is a bulk factor (each I/P converter has a volume more or less equal to 500 cm 3 ) but commercialized smaller I/P converters could be expected in the future. The supply requirement is essentially pneumatic; the electrical supply requirement is limited to the supply of the I/P converters, sensors and a robot controller. The robot pneumatic supply of the arm is ensured by the University's air-supply system (6 bar). This supply could be attached to the robot under the form of a high-pressure gas tank, although the energy range of the robot is limited by the high air consump- [START_REF] Davis | Enhanced modeling and performance in braided pneumatic muscle actuators[END_REF] or by a reprocessing of the pressurized air between the two actuator muscles.

Experimental Validation of the Prototype

The robot was tested in a teleoperation mode by means of two three-axis joysticks that control joints 1, 2, 4 and joints 5, 6, 7, and a swing bar controlling joint 3. Joint position can be controlled in open loop since the antagonistic muscle actuator is stable in open loop or in closed loop through a simple linear PID controller; the problem of actuator control is not discussed in this paper, as it was dealt with in our extensive paper [START_REF] Tondu | Modeling and control of McKibben artificial muscle robot actuators[END_REF] and discussed elsewhere (Hesselroth et al. 1994; Van der Smagt, Groen, and Schulten 1996;

Cai and Yamaura 1997). The teleoperation of the robot has a double goal: testing real joint ranges especially against gravity, which is a fundamental character for humanoid robotics, and testing effective robot compliance when tasks involving a contact of the robot with its environment are performed.

Test of Robot Joint Ranges Against Gravity

It is well known that the joints of a robot-manipulator are influenced by dynamic effects from other joints and from the static effect of gravity. In classical industrial robots, these effects are generally diminished by using high ratio speed reducers associated with high-speed rotating actuators. Such a strategy is impossible in the case of our robot because the natural compliance of the artificial muscle will be lost. Direct driving or driving by angle multiplier of moderate ratio are the only adapted driving modes. The first consequence is the dynamic non-linear behavior of the artificial muscle actuator, which must be taken into account by the robot [START_REF] Tondu | Modeling and control of McKibben artificial muscle robot actuators[END_REF]. It must also be noted, however, that the field of applications to be performed by this kind of soft robot-arm necessitates generally moderate speeds and accelerations. The second consequence that we wish to study here is the static behavior of robot against gravity. This behavior is submitted to the dependence of the actuator and gravity torques on actuator angle, as Figure 11 illustrates. Using the actuator static model of eq. ( 10), including a speed increaser of ratio m V , and notations of Figure 11(a), the following equation expressing the total joint torque submitted to the gravity effect is obtained:

T joint = (1/m V )k 1 P -(1/m V ) 2 k 2 θ joint -mgl cos θ joint . (14)
A typical simulation of joint eq. ( 14) in a [-90 • ,+90 • ] angular range is given in Figure 11(b) assuming a P constant equal to the maximum value (P 0 /2). It highlights the fact that the total torque can take a minimum value for a joint angle between 0 • and 90 • . If this minimum value is negative, due, for example, to a load at the robot end-effector, the maximum corresponding joint angle will be equal to the joint value corresponding to a total torque equal to zero. The flexion-extension and abduction-adduction shoulder joints are particularly sensitive to this situation because their actuator could have to face maximum resistive gravity torque when the forearm is in full extension. When a large joint range is desired, associated with the need for high torques, as required by the shoulder joints, it is possible to associate power muscles with a speed increaser of moderate ratio. Figure 11(c) gives the simulation result obtained with the joint 1 parameters of the 7R robot moving in full forearm extension. Without any load, it appears clearly that the full 180 • desired joint range is obtained. This is also the case if the robot tool carries a 1 kg load. However, if the end-effector robot load becomes equal to 2 kg, the effective joint range is limited to about 90 • due to the typical joint response against gravity. Such simulation illustrates a typical property of robot limbs actuated with artificial muscles: joint ranges are highly gravity-dependent.

Figure 12 shows the variations in movement of the seven joints of the robot controlled at very slow speed to eliminate dynamic effects. Joint 1 and joint 2 motions were recorded at full-arm extension. The theoretical estimated full angular range was checked for every joint and the effect of carrying loads tested. By comparison with anatomical values, all joint ranges can be performed by our anthropomorphic robot-arm except for joint 2, for which current range at full forearm extension is limited to about [0˚, 100 • ]. Such a phenomenon highlights the difficulty of designing compact and powerful shoulders performing both flexion-extension and abductionadduction. For other mechanical reasons, Honda humanoid P3 is unable to lift its extended arms higher than the horizontal and robot Asimo limits its maximum shoulder "joint mount-ing" to about 110 • , due to the overlapping of the two shoulder freedoms of movement, as mentioned in the available technical data (http://world.honda.com/ASIMO/). The paradoxical difficulty of lifting a 6 kg anthropomorphic robot-arm actuated with artificial shoulder artificial muscles generating a maximum force close to 5000 N highlights a double specificity of human body mobility. First, human body links are relatively light since the upper limb weighs only about 3.6 kg [START_REF] Winter | Anthropometry (Chapter 3) and Anthropometry, Kinematic and Force Plate Data (Appendix A). Biomechanics of Human Movement[END_REF]. Only the use of composite materials would result in our robot-arm prototype being lightened to this anatomical value. Secondly, the control of natural motions is based on muscular synergies. Let us recall that the human shoulder is governed by about 10 main muscles, for which the control sequence depends on the arm elevation angle.

Soft Task Test

To test its ability to perform soft grasping manipulation tasks, we teleoperated the robot-arm in joint space conditions to grasp a glass and put it back on to a hard surface, as Figure 13 illustrates. When the glass is on the table it is still possible, even in joint space control condition, to make it slide on the table without breaking it. Such an academic task stresses the surprising skills of artificial muscle actuators to imitate the human touch.

The compliance of the robot-tip is not controlled in this experiment but it is important to note that it is possible to define a stiffness matrix in the Cartesian space. If we note q as the vector of the robot joint variables, T r the vector of joint restoring torques, x the robot-tip location vector and F r the corresponding restoring force-torque vector in the Cartesian space, the application of the virtual works theorem gives δx T .F r = δq T .T r .

By considering the Jacobian matrix of the robot J and its pseudo-inverse J + , we could write F r = -Kδx with K = (J + ) T K q J + (15)

where K q represents the joint stiffness matrix and K is the corresponding Cartesian stiffness matrix. K q is a diagonal matrix, the terms of which have, from eq. ( 9), the general form K 2 (P 1 + P 2 ) reducing to k 2 in case of symmetrical pressure control, if necessary multiplied by a transmission term such as (1/m V ) 2 . In the framework of our paper, eq. ( 15) expresses a simple fact: if we motorize a robot-arm kinematically analogous to the human arm with actuators whose stiffnesses are similar to those of corresponding human joints, the compliance at the tip of the robot-arm will be similar to that of our arm. Further experiments will aim to quantify the ability of such 7R anthropomorphic robot-arms actuated with artificial muscles, to reproduce natural human touch in world space.

Conclusion

McKibben artificial muscles offer large possibilities for the actuation of human-friendly artificial limbs. The pneumatic actuator, defined by two identical antagonistic McKibben muscles driving a revolute axis, presents very interesting power-to-weight and power-to-volume ratios in comparison with other fluid or chemical artificial muscles. It also presents a true analogy with a skeletal muscle in generating a naturally damped contraction force decreasing with the muscle length. The consequence is a natural joint compliance, which generates the "human touch" of the robot. However, the soft materials (rubber and textile fibers) used in the muscle de-sign are a complicating factor leading to force-length characteristics difficult to model, and presenting a hysteresis phenomenon. The maximum contraction ratio is a real disadvantage for generating a wide joint range. Our research has aimed to test the validity of designing a robot-arm actuated by pneumatic artificial muscles, containing a human-like size forearm and arm with joint ranges approaching elementary physiological joint motion ranges in abduction-adduction, flexionextension, and external-internal rotation. Our 7R prototype motorized by antagonistic McKibben artificial muscles has validated this possibility; it fits into the development of soft arms, as produced by Bridgestone, the Japanese tire manufacturers, leaders in the field in the 1980s and 1990s. The "human touch" of our 7R anthropomorphic robot arm has been highlighted by controlling it by teleoperation. However, the natural joint compliance, a direct consequence of the springlike functioning of artificial muscles, implies a dependence of the actuator torque with joint angle, particularly important for shoulder control. It implies that the robot-arm is sensitive to gravity. In particular, the elevation or adduction of the robot-arm shows a particular point at which the joint torque becomes minimal against gravity. A negative value of this particular point limits the joint range, especially when the robot lifts a load. Although joint 1 and joint 2 muscles of our robot can develop a particularly high maximum force, shoulder elevation is limited in comparison with human shoulder possibilities. Better shoulder motion abilities could be obtained by further anthropomorphism. It is known that the human arm is controlled by more than 50 muscles in a complex synergy. The use of artificial muscles such as a robot actuator of a human-like robot would lead to the development of a true artificial musculature, controlling a model of the shoulder complex as the Washington Biorobotics Laboratory (http://brl.ee.washington.edu) anthroform project has tried to initiate [START_REF] Hannaford | The anthroform biorobotic arm: a system for the study of spinal circuits[END_REF].

Fig. 1 .Fig. 2 .

 12 Fig.1. Correspondence between roll-pitch-yaw mechanical notation and physiological joint motion notation. (a) Classical marine metaphor of roll-pitch-yaw motion of a mobile link seen as a boat sailing along a given reference Z-axis. (b) Physiological joint motions of a mobile link jointed to a given link: external-internal rotation corresponds to a roll motion, flexion-extension to a pitch motion and abduction-adduction to a yaw motion (note that the choice for a positive sense is dependent on link physiology).

  Fundamental 7R anthropomorphic structure defined in the joint physiology terminology: (a) zero-anatomical position; (b) considered zero-robot position (indicated joint ranges correspond to anatomical given by[START_REF] Kapandji | Physiology of the Joints[END_REF].

Fig. 4 .

 4 Fig. 4. Linear model of muscle static force production.

  Fig. 5. Isometric and isotonic recordings of McKibben muscle responses: (a) isometric force response for a 7R robot wrist muscle; (b) isometric force response for a 7R robot shoulder muscle; (c) and (d) isotonic response in pressure, force, and position x = (l 0l)/ l 0 for the robot shoulder muscle corresponding to (b) isometric response. muscle develops a maximum power of about 250 W; Bouisset

4.

  From the cylindrical model of the McKibben muscle, it can indeed be derived in accordance with the equations given in Section 3.2 that the muscle volume increases during contraction according to a slope dV (ε)/dε = (π r 2 0 l 0 )[3 cos 2 (α 0 )(1ε) 2 -1]/ sin 2 (α 0 ) positive descending.

  Fig. 6. Antagonistic McKibben muscle actuator: (a) functioning principle; (b) application of the equilibrium point principle to the antagonistic muscle actuator (illustrated in the case of positive θ actuator angle).

  Fig. 7. Mechanical design of the robot shoulder: (a) general diagram illustrating the independence principle for joint 1 and joint 2 motions; (b) photograph of the two muscle pairs driving the robot shoulder joints.

  Mechanical design of the robot wrist: (a) photograph of the robot wrist showing joint 6 direct drive; (b) diagram of joint 7 transmission by swing-bar; (c) photograph of joint 7 transmission.

Fig. 9 .

 9 Fig. 9. General diagram of the main parts of the anthropomorphic robot-arm mechanical design.

Fig. 11 .

 11 Fig. 11. Theoretical analysis of the behavior of the antagonistic artificial muscle actuator behavior against gravity: (a) notations; (b) dependence of the motor torque and the gravity torque on the joint angle; (c) simulation of the joint torque according to the load carried by the robot tool.

Fig. 12 .

 12 Fig. 12. Experimental recording of the 7R anthropomorphic joint evolutions (joint zero values correspond to the Figure 3(b) zero-robot position): (a) joint 1-3 motions performed with the arm fully extended; (b) joint 4-7 motions.

Fig. 13 .

 13 Fig. 13. Picking up and placing: the teleoperated task of a fragile object, highlighting both actuator sensitivity and its natural compliance.

Table 1 . Choice of Positive Sense of the Physiological-Type Arm Motions

 1

	Rotation

Table 2 . Actuators Characteristics

 2

					Pulley-Chain	Transmission
			Muscle		System	System			Estimated
		Characteristics	Characteristics	Characteristics		Estimated	Peak
	Joint	l 0	r 0	α 0		R			Joint Range	Torque
	Number (mm) (mm) (deg)	ε 0	(mm)	Mode	m v	(deg)	(N.m)
	1	195	12	17	0.11	20.5	Cylindrical gear	1.5	[0 • -180 • ]	61
	2	230	12	17	0.15	23	Double bevel	1.5	[0 • -100 • ]	67
							gear in
							opposition
	3	140	8.5	20	0.15	18	Bevel gear	1.5	[-100 • ,+100 • ]	18
	4	195	8.5	19	0.10	14.5	Direct drive	1	[0 • ,+150 • ]	25
	5	120	7	23	0.08	12.5	Bevel gear	2	[-90 • ,+90 • ]	5
	6	120	7	23	0.07	12.5	Direct drive	1	[-35 • ,+35]	10
	7	120	7	23	0.11	17.5	Direct drive by	1	[-45 • ,+45 • ]	13
							cross-bar
	tion of McKibben muscles, as mentioned earlier. However, the
	latter could be reduced by a pre-filling of the artificial muscle

This result is consistent with the fact that the McKibben muscle cannot contract without increasing its volume[START_REF] Tondu | Theory of an artificial pneumatic muscle and application to the modelling of McKibben artificial muscle[END_REF].

Acknowledgments

We would like to express our sincere thanks to the anonymous reviewers for their very pertinent remarks.


Samy Abbes
email: samy.abbes@univ-paris-diderot.fr

Toward Uniform Random Generation in 1-safe Petri Nets

We study the notion of uniform measure on the space of infinite executions of a 1-safe Petri net. Here, executions of 1-safe Petri nets are understood up to commutation of concurrent transitions, which introduces a challenge compared to usual transition systems. We obtain that the random generation of infinite executions reduces to the simulation of a finite state Markov chain. Algorithmic issues are discussed.

1-Introduction

Petri nets are formal models designed to describe and analyze the behavior of concurrent systems. Among the many kinds of systems where Petri nets may be introduced to formally describe a concurrent dynamics, distributed databases [9] and telecommunication networks [START_REF] Benveniste | Diagnosis of asynchronous discrete-events systems: a net unfolding approach[END_REF] are two typical examples. Both examples involve temporal evolution on the one hand, and the paradigm of resource sharing on the other hand, where resources are "spatially" distributed. Requests for resources are local, in such a way that any two actions requiring disjoint sets of resources may be considered as parallel.

Since their initial introduction in the 1960's, several variants of Petri nets have been studied. In this paper, we shall limit ourselves to 1-safe Petri nets, which we briefly define now. An unmarked Petri net is a triple N = (P, T, F ), where P and T are two finite and disjoint non empty sets of places and of transitions respectively, and F ⊆ (P × T ) ∪ (T × P ) is called the flow relation. Graphically, places are traditionally represented by circles and transitions are represented by squares or rectangles (see Figure 1). The flow relation is depicted by arrows from places to transitions and from transitions to places.

Given a transition t ∈ T , the preset • t and the postset t • of t are the sets of places defined as follows:

• t = {p ∈ P : (p, t) ∈ F } , t • = {p ∈ P : (t, p) ∈ F } .

It is assumed that the preset and the postset of any transition are both non empty.

A marking of N is any integer valued function M : P → N. The marking is said to be 1-safe, or simply safe, whenever M (p) ≤ 1 for all places p ∈ P . The number M (p) is interpreted as a number of tokens lying in the place p. Tokens are graphically figured inside places, as in Figure 1. Given a marking M of N , and a transition t ∈ T , we say that t can fire from M , or that M enables t, whenever M (•) > 0 on • t .

If t can fire from M , then the firing rule M t -→ M defines the new marking M as follows (see the commentary below and the illustration in Figure 1):

∀p ∈ P M (p) =          M (p), if p / ∈ ( • t ∪ t • ) M (p) -1, if p ∈ • t \ t • M (p) + 1, if p ∈ t • \ • t M (p), if p ∈ ( • t ∩ t • ) $ $ • Ö Ö $ $ • Ö Ö • a b firing of transition b ------------------→ $ $ • Ö Ö $ $ Ö Ö • a b • Figure 1: Illustrating the firing rule in Petri nets
The intuitive interpretation of the firing rule is as follows. The tokens located in places of the preset • t are resources, needed for firing and consumed by the transition when fired; this explains the second rule. The firing of t also produces new resources, i.e., tokens; the third rule specifies that the new tokens are created in the postset of t. If a place belongs to both the pre-and the postset of t, then we can see the fourth rule as the simultaneous instance of both the second and the third rule: an existing resource is consumed and immediately after a new one is created, replacing the first one. Finally, the first rule specifies that only places "near" the firing transition are concerned by the firing rule. Note that the total number of tokens in the net may differ once the transition has fired.

Given a marking M 0 , a firing sequence from M 0 is a sequence (t 1 , . . . , t k ) of transitions such that, for some markings M 1 , . . . , M k defined inductively, the firing rules M i-1 ti -→ M i hold for all i = 1, . . . , k . Of course, if such markings exist, they are unique. Retaining only the last marking, we introduce the obvious notation

M 0 t1•••t k ----→ M k , with the convention M 0 ε -→ M 0
for the empty sequence ε, which is considered as a firing sequence.

Any marking M such that M 0

t1•••t k ----→ M holds for some firing sequence t 1 • • • t k is said to be reachable from M 0 .
A Petri net is a quadruple N = (P, T, F, M 0 ), such that (P, T, F ) is an unmarked Petri net, of which M 0 is a marking, called the initial marking of the net. The Petri net N is said to be 1-safe, or simply safe, whenever all markings reachable from M 0 are safe; this implies of course that M 0 itself is safe. By convention, we assume that M 0 is fixed once and for all, and by a marking M of the net we mean any marking reachable from M 0 .

Based on firing sequences, authors have introduced a refined point of view on the executions of Petri nets. The idea is to underline the concurrency features of the model. Indeed, the firing rule enlightens that a special status should be given to pairs (t, t ) of transitions such that (

• t ∪ t • ) ∩ ( • t ∪ t • ) = ∅;
in this paper, we shall say that two such transitions are distant. According to the firing rule, if both t and t are enabled by some marking M , and if they are distant, then: 1) the firing of one does not prevent the firing of the other one, and 2) the markings resulting from the two firing sequences t • t and t • t, fired from M , are the same.

Henceforth, instead of mere firing sequences, it is natural to consider equivalence classes of firing sequences, with respect to some congruence R, such that (t • t , t • t) ∈ R for any two distant transitions t and t . We define thus R as the smallest congruence on the set firing sequences that contains all pairs (t • t , t • t), for t and t two distant transitions. Hence, two firing sequences are congruent whenever one can pass from one to another by applying a finite number of times elementary transformations of the form:

(x) • t • t • (y) - → (x) • t • t • (y)
where (x) and (y) are any two words on the alphabet T , and t and t are two distant transitions.

Equivalence classes of firing sequences are called configurations; intuitively, they only capture the causal relations between transitions, and they leave unspecified the chronological relations between distant and unrelated transitions; see for instance G. Winskel's notion of event structure and of unfolding of Petri net for a precise account on this interpretation [START_REF] Nielsen | Petri nets, event structures and domains, part 1[END_REF].

In this paper, we are interested in the random sampling of configurations of a given safe Petri net. We insist that it shall not be confused with the random sampling of firing sequences. The later ultimately reduces to the standard procedure of sampling in a large, but finite transition system, with the set of markings of the net as set of states. By contrast, the sampling of configurations is by no means standard-at least, not on first sight.

Therefore, we consider configurations as our primary objects of interest, instead of their underlying firing sequences. It is consistent with the general idea that the concurrent firing of distant transitions should not be ordered, since the ordering of parallel transitions does not correspond to any "observable" of the system and thus would be rather artificial.

In order to properly deal with configurations, we rely first on the notion of trace monoid (also called heap monoid [START_REF] Viennot | Heaps of pieces, I : basic definitions and combinatorial lemmas[END_REF] or free partially commutative monoid [START_REF] Cartier | Problèmes combinatoires de commutation et réarrangements[END_REF]). More precisely, the exact notion that we need is that of a trace monoid acting on a finite set. Associated with the action of a trace monoid, is a notion of uniform probability measure on the space of "infinite traces", developed in a previous work [START_REF] Abbes | Dynamics of asynchronous systems[END_REF]. In the context of Petri nets, infinite traces correspond to the infinite executions of the net.

The purpose of this paper is thus twofold: 1) give an account on the uniform measure relative to the action of a trace monoid on a finite set; and 2) analyze the application of this theory in the context of safe Petri nets and for random sampling purposes. The intended applications related to random generation are of two kinds: random generation of infinite executions of a Petri net, uniformly distributed and once this has been properly defined, on the one hand; random generation of a finite execution, uniformly distributed among all finite executions of size n, on the other hand. The results of the paper allow to tackle the first task, at least theoretically, although some algorithmic issues inevitably remain and are discussed. The second task, concerning finite sampling, is not discussed here because of space constraints. Let us mention however that the solution to a similar, but simpler problem treated in [START_REF] Abbes | Uniform generation in trace monoids[END_REF], could probably be extended to this case with few modifications. In turn, these are the basic building blocks needed, for instance, for probabilistic model checking of formal properties of concurrent systems.

The application to Petri nets will be treated on an example. Henceforth, the contributions of this paper consist mainly in the new questions that we formulate concerning the uniform measure relatively to the specific model of Petri nets, compared to more general asynchronous systems.

We underline that the probabilistic dynamics to which we are naturally brought radically differs from that of stochastic Petri nets [START_REF] Baccelli | Synchronization and Linearity[END_REF]-this is due to the radical turn we have taken in our analysis by considering configurations instead of firing sequences, as explained above.

The outline of the paper is as follows. Section 2 recalls the basics on trace monoids and their combinatorics. Section 3 is devoted to the action of trace monoids and their combinatorics, and introduces the associated notion of uniform measure. It is only in Section 4 that we meet with Petri nets again, applying the results from Section 3 to their case.

2-Trace monoids and their boundary

Let Σ be an alphabet, that is to say, a finite set, the elements of which are called letters. Let also I ⊂ Σ × Σ be a binary relation, that we assume to be symmetric and irreflexive, and that is called an independence relation. Let Σ * denote the monoid of Σ-words, and let R be the smallest congruence on Σ * containing all pairs (ab, ba) for (a, b) ∈ I. The trace monoid M = M(Σ, I) is the quotient monoid M = Σ * /R . In other words, M is the presented monoid:

M = Σ | ab = ba for all (a, b) ∈ I .
Elements of M are called traces. The concatenation of traces is denoted with the dot "•", and the identity element is denoted ε. Clearly, if x and y are two R-congruent words, they have the same length. This defines a length function | • | : M → N on traces. Clearly also, the left divisibility relation on M, defined by x ≤ y ⇐⇒ ∃z y = x • z, is a partial order on M.

A clique of M is any product of the form γ = a 1 • . . . • a k , where a 1 , . . . , a k are letters such that i = j =⇒ (a i , a j ) ∈ I . Note in particular that the a i 's are pairwise distinct. Let C denote the set of cliques; this is a finite set. Let also C = C \ {ε} denote the set of non empty cliques. A pair (γ, γ ) ∈ C × C is said to be in normal form whenever, for each letter b occurring in γ , there is a letter a occurring in γ such that (a, b) / ∈ I. We denote it by the symbol γ → γ . A sequence (γ 1 , . . . , γ n ) in C is said to be normal whenever γ i → γ i+1 holds for all i ∈ {1, . . . , n -1} . It is well known that traces admit the following normal form [START_REF] Cartier | Problèmes combinatoires de commutation et réarrangements[END_REF][START_REF] Viennot | Heaps of pieces, I : basic definitions and combinatorial lemmas[END_REF]: for every x ∈ M \ {ε}, there exists a unique integer n ≥ 1 and a unique normal sequence (γ 1 , . . . , γ n ) in

C such that x = γ 1 • . . . • γ n .
The existence of this normal form is the basis for establishing the following combinatorial results. Let µ(z) be the polynomial, called Möbius polynomial, and defined by:

µ(z) = γ∈C (-1) |γ| z |γ| .
Let also the growth series G(z) be defined by:

G(z) = k≥0 λ(k)z k , λ(k) = #S(k) , S(k) = {x ∈ M : |x| = k} .
Then G(z) is a rational series, inverse of the Möbius polynomial: G(z) = 1/µ(z) . Furthermore [START_REF] Goldwurm | Clique polynomials have a unique root of smallest modulus[END_REF][START_REF] Krob | Computing the average parallelism in trace monoids[END_REF][START_REF] Csikvári | Note on the smallest root of the independence polynomial[END_REF], µ(z) has a unique root of smallest modulus. This root, say p 0 , is real and lies in (0, 1], and coincides with the radius of convergence of the series G(z). We shall extend in next section this well known result in the framework of a trace monoid acting on a finite set.

The existence of the normal form entails that traces are in bijection with finite paths in the finite graph of non empty cliques (C, →). It is thus natural to define infinite traces as infinite paths in the very same graph (C, →). Let ∂M denote the set of infinite traces. The set ∂M is called the boundary of M. It is standard that the set ∂M, equipped with the natural topology, is metrisable and compact.

It is a bit less standard to extend the partial ordering relation ≤ from M to ∂M. Let ξ = (γ 1 , γ 2 , . . .) be an infinite trace, and let x ∈ M. Then we put x ≤ ξ if and only if x ≤ (γ 1 • . . . • γ k ) in M for all k large enough. And we define the visual cylinder ↑ x, as the following subset of ∂M :

↑ x = {ξ ∈ ∂M : x ≤ ξ} .
We will be interested in probability measures on the space ∂M equipped with its Borel σ-algebra. The following result shows the importance of visual cylinders in this respect (see [2, § 2.2]).

• Lemma 2.1-Any probability measure m on ∂M is entirely characterized by the countable collection m( ↑ x), for x ranging over M.

We conclude this section by illustrating the above notions on a concrete example. It will be given a Petri net interpretation below in Section 4.

Let Σ = {a, b, c, d, e} and let M = Σ | ad = da, ae = ea, bd = db, be = eb, ce = ce . The set of non empty cliques is C = {a, b, c, d, e, ad, ae, bd, be, ce} . For instance, the trace adebc = daebc = dabec has the normal form (ad) → (be) → (c).

The Möbius polynomial of M is: µ(z) = 1 -5z + 5z 2 , and its root of smallest modulus is

p 0 = 1 2 -1 2 √ 5 .

3-Uniform measure relative to the action of a trace monoid

Let M = M(Σ, I) be a trace monoid. Let also Y be a finite set of states. We assume given a right monoid action of M on Y , that is to say, a function ϕ : Y × M → M denoted by ϕ(s, x) = s • x, satisfying the following:

∀s ∈ Y s • ε = s , ∀(s, x, y) ∈ Y × M × M s • (x • y) = (s • x) • y .
The pair (Y, M), with the action understood, is called an asynchronous system. The traces x ∈ M are thought of as actions featuring parallelism; the effect of the action x ∈ M on a state s ∈ Y is to change the state of the system into the new state s • x .

A desirable feature is to disable some actions, depending on the current state of the system. We render this feature by assuming the existence of a particular state ⊥, such that ⊥ • x = ⊥ for all x ∈ M. We put Y = X ∪ {⊥} with ⊥ / ∈ X, so that the "real" states are actually the elements of X; we will always restrict our attention to those pairs (s, x) ∈ X × M such that s • x = ⊥. Hence, it is customary to introduce the following notation:

∀s ∈ X M s = {x ∈ M : s • x = ⊥} .
We say that the action is irreducible if, for every pair (s, t) ∈ X × X of states, there exists an action x ∈ M such that x = ε and t = s • x. We shall always assume that the action under consideration is irreducible.

On the combinatorics side, the analogous of the growth series is the following matrix of formal series:

G = (G s,t ) (s,t)∈X×X , G s,t (z) = x∈M : x•s=t z |x| .
Note that each line of G, say indiced by s ∈ S, sums up to the growth series G s of the subset M s , defined by:

G s (z) = k≥0 λ s (k) , λ s (k) = #{x ∈ M s : |x| = k} . (1)
Let us define the Möbius matrix M by:

M = (M s,t ) (s,t)∈X×X , M s,t (z) = γ∈C : s•γ=t (-1) |γ| z γ ,
where C denotes the set of cliques of M . Then M (z) is the formal inverse of the growth matrix: G(z)M (z) = I, where I is the identity matrix of size |X|× |X| (see [START_REF] Abbes | Dynamics of asynchronous systems[END_REF]Theorem 5.10]). Given our hypothesis that the action is irreducible, all entries in the growth matrix G have the same radius of convergence, say q 0 ∈ (0, 1]. It is also the common radius of convergence of the growth series G s (z), for s ∈ X, defined in (1). We call q 0 the characteristic root of the system (Y, M).

Let the theta polynomial be the polynomial with integer coefficients defined by θ(z) = det M (z). Then q 0 is the smallest positive root of θ(z) [START_REF] Abbes | Dynamics of asynchronous systems[END_REF]Theorem 5.11].

The above elements of combinatorics allow us to introduce a notion of uniform measure, as follows. For each state s ∈ X, and for each real number q ∈ (0, q 0 ), let m s,q be the discrete probability measure on M defined by:

m q,s = 1 G s (q) x∈Ms q |x| δ {x} ,
where δ {x} is the Dirac measure at x. Then, according to [1, Theorem 5.16], the family (m s,q ) q<q0 converges weakly, as q → q - 0 , toward a probability measure ν s concentrated on ∂M, and which satisfies the following property (recall that a probability measure on ∂M is entirely determined by its values on visual cylinders):

∀s ∈ X ∀x ∈ M s ν s ( ↑ x) = q |x| 0 Γ(s, s • x) , (2)
and ν s ( ↑ x) = 0 whenever x / ∈ M s , i.e., whenever s • x = ⊥ . In (2), Γ(•, •) : X × X → (0, ∞) is a function satisfying the following cocycle identity: Γ(s, u) = Γ(s, t)Γ(t, u) for all s, t, u ∈ X .

We call the family ν = (ν s ) s∈X the uniform measure of the asynchronous system (Y, M), although it is only each individual ν s which is actually a probability measure. Clearly, the state s indexing ν s is to be seen as the initial state of the system. The form (2), together with the cocycle identity, immediately yields the following property:

∀(s, x, y) ∈ X × M × M ν s ↑ (x • y) = ν s (x) • ν s•x ( ↑ y) . (3)
This chain rule (3) justifies that the uniform measure ν is called a Markov measure; but there are other families satisfying (3) than the uniform measure that we just constructed.

Regarding a uniqueness result, we leave the following question open: is there a unique pair (q 0 , Γ), where q 0 ∈ (0, 1] and Γ(•, •) : X × X → (0, ∞) satisfies the cocycle identity, such that the formula ν s ( ↑ x) = q |x| 0 Γ(x, s • x) for x ∈ M s and ν s ( ↑ x) = 0 for x / ∈ M s , defines a probability measure on ∂M for all s ∈ X?

From the sampling point of view, the following realization result is interesting. Let s 0 ∈ X be a fixed initial state, and let ξ ∈ ∂M be an infinite trace distributed according to ν s0 . Let ξ be given as the infinite path ξ = (C 1 , C 2 , . . .) in (C, →) (see Section 2). Here, we see (C k ) k≥1 as a sequence of random variables under the probability ν s0 . Finally, for k ≥ 0, let S k be the state of the system reached after the action of the k first cliques:

S k = s 0 • (C 1 • . . . • C k ) , with S 0 = s 0 by convention.
Then, according to [1, Theorem 4.5], under ν s0 , the sequence of pairs (S k-1 , C k ) k≥1 is a homogeneous Markov chain with values in (X × C), that we call the Markov chain of states-and-cliques. Furthermore: 1) the transition matrix of the chain is independent of s 0 ; and 2) both the initial measure and the transition matrix of the chain have algebraic expressions involving only q 0 and Γ(•, •). See [START_REF] Abbes | Dynamics of asynchronous systems[END_REF] for the explicit expressions.

Hence, both for theoretical and for sampling purposes, the problem comes down to determining the characteristic root q 0 and the cocycle Γ(•, •). As for q 0 , its characterization as a particular root of the theta polynomial yields the following algorithmic issues: 1) list all cliques of M, which is very well known to be hard (exponential in the size of Σ in worst case [START_REF] Tomita | The worst-case time complexity for generating all maximal cliques and computational experiments[END_REF]); 2) compute a determinant of size |X| × |X|, the entries of which are polynomials of degree at most α, where α is the maximal size of a clique; 3) solve the equation θ(z) = 0 for the smallest positive root, where θ is of degree at most α|X| and α is the maximal size of a clique.

Regarding the cocycle Γ(•, •) , we have not found a systematic way of determining it. We will see in next section how to determine it in practice on an example. Based on the uniqueness conjecture stated above, we can expect that this technique actually works in general, yielding an algorithmic way of obtaining the cocycle.

4-Application to 1-safe Petri nets

In this section, we study 1-safe Petri nets as asynchronous systems as defined in the previous section. Although this point of view is quite natural, it does not seem to have been adopted by other authors yet, not formally at least.

Let N = (P, T, F, M 0 ) be a 1-safe Petri nets. Define the set of states X as the state of reachable markings of N . We consider Σ = T as alphabet, and as independence relation I, we take

I = {(t, t ) ∈ T × T : ( • t ∪ t • ) ∪ ( • t ∪ t • ) = ∅},
hence the set of distant pairs as defined in Section 1. Finally, let M = M(Σ, I) be the associated trace monoid (see Section 2).

Putting Y = X ∪ {⊥}, we define an action of M on Y as follows. Let M ∈ X, and let x ∈ M, of which (t 1 , . . . , t k ) is a representative sequence of transitions. Then, if (t 1 , . . . , t k ) is a firing sequence from M , we put M • x = M such that M t1•...•t k -----→ M ; and if (t 1 , . . . , t k ) is not a firing sequence from M , we put M • x = ⊥. Obviously, this definition does not depend on the representative sequence (t 1 , . . . , t k ). It is clear also that M • (x • y) = (M • x) • y holds for all x, y ∈ M.

Let us consider for instance the Petri net depicted in Figure 2, with the initial marking, say M 0 , depicted. Clearly, this net is safe, and its only reachable marking, say M 1 , other than M 0 is the marking obtained from M 0 by firing transition b. The associated trace monoid is M = a, b, c, d, e, | ad = da, ae = ea, bd = db, be = eb, ce = ec , that is to say, the example trace monoid introduced at the end of Section 2. The action obeys the following description, which implies in particular that it is an irreducible action:

M 0 • a = M 0 M 0 • b = M 1 M 0 • c = ⊥ M 0 • d = M 0 M 0 • e = M 0 M 1 • a = ⊥ M 1 • b = ⊥ M 1 • c = M 0 M 1 • d = M 1 M 1 • e = M 1
Note however that this description does not characterize the action by itself, since it does not render the concurrency features encoded in the independence relation I.

The Möbius matrix and the theta polynomials are:

M (z) = 1 -3z + 2z 2 -z + 2z 2 -z + z 2 1 -2z θ(z) = (1 -z)(1 -2z)(1 -2z -z 2 ) .
The characteristic root is thus Let ν = (ν M0 , ν M1 ) be the uniform measure associated with the action. Let also Γ(•, •) be the cocycle associated with ν, so that ν M ( ↑ x) = q |x| 0 Γ(M, M • x) for all M ∈ {M 0 , M 1 } and for all x ∈ M M . Then, by the cocycle property of Γ, one has Γ(M 0 , M 0 ) = Γ(M 1 , M 1 ) = 1 on the one hand, and one has Γ(M 1 , M 0 ) = 1/Γ(M 0 , M 1 ) on the other hand. Hence Γ(•, •) is entirely determined by the single value λ = Γ(M 0 , M 1 ).

q 0 = √ 2 -1. Ò Ò ( ( • R R a b 0 0 Ò Ò ( ( • Ò Ò ( ( • v v R R c v v u u d e
In order to retrieve the value λ, we write down the following equation. This is equivalent [START_REF] Abbes | Uniform and Bernoulli measures on the boundary of trace monoids[END_REF][START_REF] Abbes | Dynamics of asynchronous systems[END_REF] to writing down the total probability law for the first clique C 1 , starting from marking M 0 :

γ∈C : M0•γ =⊥ (-1) γ ν M0 ( ↑ γ) = 0 .
Since ν M0 ( ↑ γ) = q |γ| 0 Γ(M 0 , M 0 • γ) on the one hand, and since the cliques γ such that M • γ = ⊥ range over {ε, a, b, d, e, ad, ae, bd, be}, this writes as:

1 -3q 0 Γ(M 0 , M 0 ) -q 0 Γ(M 0 , M 1 ) + 2q 2 0 Γ(M 0 , M 0 ) + 2q 2 0 Γ(M 0 , M 1 ) = 0 .
Since Γ(M 0 , M 0 ) = 1, introducing the unknown λ = Γ(M 0 , M 1 ) yields:

λ = 1 -3q 0 + 2q 2 0 q 0 (1 -2q 0 ) = -1 + 3q 0 1 -2q 0 ,
the later equality following from the equation 1 -2q 0 -q 2 0 = 0. Introducing the value q 0 = √ 2 -1 finally yields the simple expression: λ = √ 2. Therefore the visual cylinders are given probabilities with the following simple expressions:

∀x ∈ M M0 ν M0 ( ↑ x) = ( √ 2 -1) |x| , if M 0 • x = M 0 √ 2( √ 2 -1) |x| , if M 0 • x = M 1 ∀x ∈ M M1 ν M1 ( ↑ x) = ( √ 2 -1) |x| , if M 1 • x = M 1 1 √ 2 ( √ 2 -1) |x| , if M 1 • x = M 0
The above formulas however are not operational for sampling purposes. One way to obtain a sampling method is to consider the Markov chain of states-and-cliques, that is to say, the random sequence of pairs (M k-1 , C k ) k≥1 , where (C k ) k≥1 is the sequence of cliques forming an infinite trace ξ and M k is the marking reached after the k th clique (see Section 3). Based on the formulas of [START_REF] Abbes | Dynamics of asynchronous systems[END_REF] and on the values found above for q 0 and for Γ(•, •), we find that the transition matrix of (M k-1 , C k ) k≥1 is the following stochastic matrix-keeping in mind that only the reachable pairs (M, γ) ∈ X × C have to be considered, see the remark below:

             -1 + √ 2 0 0 0 0 2 - √ 2 0 0 0 0 -2 + 3 2 √ 2 0 1 -1 2 √ 2 0 0 3 -2 √ 2 0 0 -1 + √ 2 0 -7 + 5 √ 2 0 3 -2 √ 2 3 -2 √ 2 0 10 -7 √ 2 0 0 -4 + 3 √ 2 -4 + 3 √ 2 -7 + 5 √ 2 0 3 -2 √ 2 3 -2 √ 2 0 10 -7 √ 2 0 0 -4 + 3 √ 2 -4 + 3 √ 2 -7 + 5 √ 2 0 3 -2 √ 2 3 -2 √ 2 0 10 -7 √ 2 0 0 -4 + 3 √ 2 -4 + 3 √ 2 0 3 -2 √ 2 0 0 -2 + 3 2 √ 2 0 -1 + √ 2 1 -1 2 √ 2 0 0 0 3 -2 √ 2 0 0 -2 + 3 2 √ 2 0 -1 + √ 2 1 -1 2 √ 2 0 0 0 0 0 0 0 0 2 - √ 2 -1 + √ 2 0 0 0 3 -2 √ 2 0 0 -2 + 3 2 √ 2 0 -1 + √ 2 1 -1 2 √ 2 0 0 0 3 -2 √ 2 0 0 -2 + 3 2 √ 2 0 -1 + √ 2 1 -1 2 √ 2 0 0             
Furthermore, the initial law of the chain, i.e., the law of (M 0 , C 1 ) , has the form δ {M0} ⊗ κ , where κ is the probability distribution on C given by:

κ(a) = -7 + 5 √ 2 κ(b) = 10 -7 √ 2 κ(c) = 0 κ(d) = 0 κ(e) = 0 κ(ad) = 3 -2 √ 2 κ(ae) = 3 -2 √ 2 κ(bd) = -4 + 3 √ 2 κ(be) = -4 + 3 √ 2 κ(ce) = 0
A remark on these values: the equality κ(c) = 0 is normal, since c is not enabled at M 0 . On the contrary, both d and e are enabled at M 0 , and yet we find that κ vanishes on d and e. This is obtained from a straightforward computation based on the formulas of [START_REF] Abbes | Dynamics of asynchronous systems[END_REF] applied to our case. But it is also worth mentioning an intuitive explanation, as follows. Assume for instance that C 1 = d . Then the system stays in M 0 , and since c is not enabled at M 0 , the next clique will be either d or e, but will not contain either a nor b, by the definition of the normal form. The same will happen for the next clique, and so on. Henceforth a and b will never occur; but this has zero probability, and it explains that κ(d) = 0. The same holds for κ(e) = 0. The same reasoning also explains that the pairs (M 0 , d) and (M 0 , e) are not reachable states of the Markov chain of states-and-cliques, and thus do not appear as lines in the transition matrix above.

Given the initial distribution and the transition matrix, one is of course able to simulate the Markov chain (M k-1 , C k ) k≥1 . In turn, the random trace (C 1 • . . . • C k ) corresponds to the k first cliques of a uniformly distributed infinite execution of the net.

We observe that the transition matrix displayed above has the following property: for any clique γ, the law of the next marking starting from (M 0 , γ) is independent of γ; and similarly when starting from (M 1 , γ). We deduce that the sequence (M k ) k≥0 is itself a Markov chain, with transition matrix:

-1 + √ 2 2 - √ 2 1 -1 2 √ 2 1 2 √ 2
Based on the results obtained for the above specific example, the following questions arise naturally regarding how to deal with a general 1-safe Petri net.

1. What would be an algorithmic way of computing the cocycle Γ(•, •)? 2. The uniform measure has the simple characterization ν s ( ↑ x) = q |x| 0 Γ(s, s • x) . By contrast, its realization as a Markov chain involves a set of states which grows, in worst case, exponentially fast with the number of places and with the number of transitions of the net. It is therefore natural to seek for other means for uniform generation than the computation of the Markov chain of states-and-cliques. This relates with the next point.

3. In our specific example, the uniform measure clearly shows a probabilistic independence between the two parallel choices between a over b on the one hand, and between d over e on the other hand. Until which extend can we generalize this spatial independence property for a safe Petri net? Does it take a more specific formulation for certain particular classes of safe Petri nets, such as freechoice nets for instance [START_REF] Best | Structure Theory of Petri nets: the Free Choice Hiatus[END_REF]? Can we use this independence property to devise an efficient sampling algorithm for executions of nets? 4. For our specific example, we have observed that the sequence of markings (M k ) k≥0 reached by the Markov chain of states-and-cliques (M k-1 , C k ) k≥1 is itself a Markov chain; this is not the case for general asynchronous systems.

Can we obtain sufficient conditions on the structure of a 1-safe Petri nets for this to hold?

Figure 2 :

 2 Figure 2: Example of a safe Petri net. The underlying trace monoid is the one introduced at the end of Section 2.