Zane Zake
email: zane.zake@ls2n.fr

Nicolò Pedemonte
email: nicolo.pedemonte@irt-jules-verne.fr

Stéphane Caro

Visual Servoing of Cable-Driven Parallel Robots with Tension Management

Cable-driven parallel robots (CDPRs) are a type of parallel robots, where cables are used instead of rigid links. This leads to many advantages, such as large workspace, low mass in motion and simple reconfiguration. The drawbacks are accuracy issues and complex cable management. Indeed, it is usual that cables become slack. That can be caused by, for example, cable mass, uncertainties in the system, and a higher number of cables than the number of degrees of freedom of the moving-platform. This reduces CDPR stiffness and degree of actuation. While visual servoing provides good accuracy and is robust to different perturbations in the system and to modeling errors, it does not deal with cable slackness. Thus, a CDPR with visual servoing can become underactuated due to cable slack. We propose in this paper to enrich visual servoing with a tension correction algorithm. Experimental results show reduction of slackness and thus avoiding slacknessrelated trajectory perturbations and loss of stability.

I. INTRODUCTION

Parallel robots have multiple kinematic chains connecting the end-effector, also known as the moving-platform (MP), and the base. In Cable-Driven Parallel Robots (CDPRs) the rigid links are replaced by cables. This change leads to many advantages, such as reduced mass in motion, large workspace, high payload capacity, and reconfigurability [START_REF] Gagliardini | Discrete Reconfiguration Planning for Cable-Driven Parallel Robots[END_REF]. The drawbacks include unidirectional cable force. Indeed, a cable can only pull, yet it cannot push. Due to this, at least m = n + 1 cables are needed to fully constrain n degrees of freedom (DoF) [START_REF] Roberts | On the inverse kinematics, statics, and fault tolerance of cable-suspended robots[END_REF], making a CDPR redundant in actuation to degree r = m -n. Generally, eight cables are used for six DoF motion, because this increases the CDPR workspace [START_REF] Gouttefarde | A Versatile Tension Distribution Algorithm for 6-DOF Parallel Robots Driven by n+2 Cables[END_REF]. In model-based approaches CDPR accuracy is directly tied to the choice of CDPR model. For instance, cable elasticity and pulley kinematics need to be taken into account to ensure acceptable accuracy [START_REF] Schmidt | Modeling Techniques and Reliable Real-Time Implementation of Kinematics for Cable-Driven Parallel Robots using Polymer Fiber Cables[END_REF][START_REF] Picard | Pulleys and Force Sensors Influence on Payload Estimation of Cable-Driven Parallel Robots[END_REF].

There are multiple ways to use a vision system for controlling CDPRs. A straightforward way is to use a camera mounted on the base to retrieve the MP pose [START_REF] Dallej | Vision-based modeling and control of large-dimension cable-driven parallel robots[END_REF]- [START_REF] Zavatta | A vision-based referencing procedure for cable-driven parallel manipulators[END_REF], which avoids the need of the Direct Kinematic Model (DKM). For CDPRs with a large workspace, multiple cameras are usually required, e.g. three cameras in [START_REF] Dallej | Modeling and Vision-Based Control of Large-Dimension Cable-Driven Parallel Robots Using a Multiple-Camera Setup[END_REF], four cameras in [START_REF] Dallej | Vision-based modeling and control of large-dimension cable-driven parallel robots[END_REF] and 44321 Nantes, France, stephane.caro@ls2n.fr six infrared cameras in [START_REF] Chellal | A Kinematic Vision-Based Position Control of a 6-DoF Cable-Driven Parallel Robot[END_REF]. When the goal is to improve the robot accuracy with respect to a target, a camera is usually mounted on the MP [START_REF] Ramadour | Grasping Objects With a Cable-Driven Parallel Robot Designed for Transfer Operation by Visual Servoing[END_REF]- [START_REF] Zake | Robust 2½D Visual Servoing of a Cable-Driven Parallel Robot Thanks to Trajectory Tracking[END_REF]. Thus, both the camera and the MP approach the target simultaneously. In that case, the MP is not observed, but its pose is required in the control scheme. For instance, the DKM was used in [START_REF] Ramadour | Grasping Objects With a Cable-Driven Parallel Robot Designed for Transfer Operation by Visual Servoing[END_REF] to retrieve the MP pose of a simplified spatial CDPR with three translational DoF. On the contrary, in [START_REF] Zake | Vision-Based Control and Stability Analysis of a Cable-Driven Parallel Robot[END_REF] and [START_REF] Zake | Robust 2½D Visual Servoing of a Cable-Driven Parallel Robot Thanks to Trajectory Tracking[END_REF] the use of DKM was avoided by coarsely estimating the MP pose. Indeed, it is possible to estimate the current MP pose, given the previous pose, the robot velocity sent to the low-level controller and the time interval, assuming that the very first initial pose is known. While visual servoing (VS) of a CDPR is very robust to many kinds of perturbations [START_REF] Zake | Vision-Based Control and Stability Analysis of a Cable-Driven Parallel Robot[END_REF][START_REF] Zake | Robust 2½D Visual Servoing of a Cable-Driven Parallel Robot Thanks to Trajectory Tracking[END_REF] and remains accurate with respect to (wrt) the target, the cables tend to become slack. Furthermore, a vicious circle can be observed. As cables become slack, they do not produce the desired MP motion. Therefore, the error in the MP pose estimation increases. The worse the MP pose estimation, the slacker the cables. Conversely, the slacker the cables, the worse the MP pose estimation. This is also true for the cable length-based DKM computation for a simplified CDPR model as in [START_REF] Ramadour | Grasping Objects With a Cable-Driven Parallel Robot Designed for Transfer Operation by Visual Servoing[END_REF].

Cable slackness is not specific to VS of CDPRs with a camera embedded on their MP. On the contrary, it is a common problem, which usually occurs over time due to discrepancies between the CDPR and its model. For example, a control technique that relies on a good knowledge of initial MP pose and then controls the robot via motor position commands is subject to cable slackness. In [START_REF] Kraus | System identification and cable force control for a cable-driven parallel robot with industrial servo drives[END_REF] the cable lengths are computed for a desired MP pose of a fullyconstrained CDPR. To avoid cable slackness and to improve robot stiffness, a complex cable length adjustment based on current tension sensor measurements and the output of a tension distribution algorithm (TDA) is implemented. The TDA can also be used in tension or torque control of a fully-constrained CDPR [START_REF] Fortin-Côté | A tension distribution algorithm for cable-driven parallel robots operating beyond their wrenchfeasible workspace[END_REF][START_REF] Picard | Stiffness Oriented Tension Distribution Algorithm for Cable-Driven Parallel Robots[END_REF], thus ensuring that tensions in cables correspond to the desired MP pose. Even with an external measurement of the MP pose, cable slackness can occur, if the model does not correspond to the real CDPR, thus often complex cable models are used [START_REF] Schmidt | Modeling Techniques and Reliable Real-Time Implementation of Kinematics for Cable-Driven Parallel Robots using Polymer Fiber Cables[END_REF][START_REF] Dallej | Vision-based modeling and control of large-dimension cable-driven parallel robots[END_REF].

The goal of this paper is to propose a simple algorithm that deals with cable slackness. Unlike the mentioned approaches, the proposed algorithm does not rely on the knowledge of the MP pose, the use of a TDA, nor a complex cable modeling depending on cable angles perceived by cameras [START_REF] Dallej | Vision-based modeling and control of large-dimension cable-driven parallel robots[END_REF].

For a redundantly actuated CDPR, having r cables slack does not prevent the MP from remaining in a given pose in static configuration. However, no matter the control ap- • F b , Fp, Fc, Fo denote the base, MP, camera and object frames respectively (resp.).

• i T j = i R j i t j 0 1
is the transformation matrix from F i to F j , where i R j and i t j are the rotation matrix and translation vector, resp.

• A -1 and A † are the inverse and the pseudo-inverse of A, resp.

• A and ê are the estimations of A and e, resp.

• s * is the desired value of s.

• [e] × denotes the cross-product matrix of vector e.

• i a is the vector a expressed in F i • A i and B i are the exit and the anchor points of the ith cable, resp. proach, cable slackness leads to the CDPR being only partly responsive to the control output and thus not behaving as expected. Moreover, along a trajectory, cable slack will be transferred to different cables during the transition between two six-cable configurations [START_REF] Merlet | Simulation of discrete-time controlled cable-driven parallel robots on a trajectory[END_REF], the MP becoming locally underactuated.

Thus, it is important to avoid cable slackness in order to improve the CDPR behavior and to keep it as stiff as possible. A tension correction algorithm is proposed in this paper to enrich VS control. Lyapunov stability analysis is used to evaluate the effect of cable slackness on the system.

The paper is organized as follows. Notations used in this paper are shown in Table I. VS of CDPRs is recalled in Section II. Tension correction is presented in Section III. Section IV contains the experimental validation and stability analysis. Finally, conclusions are drawn in Section V.

II. VISUAL SERVOING OF A CDPR

The VS control scheme of a CDPR is shown in Fig. 1. At every iteration of the control scheme an image is received from the onboard camera and it is sent to a computer vision algorithm, which outputs the current value of the feature vector s. The nature of s is a function of the chosen VS approach. For example, it can simply be the target pose expressed in the camera frame F c when a Pose-Based Visual Servoing (PBVS) is implemented [START_REF] Chaumette | Visual Servoing[END_REF]. In a more efficient VS approach, named 2½D visual servoing (2½D VS), the feature vector s contains both 2D and 3D features. For example [START_REF] Chaumette | 2 1/2 D visual servoing: a possible solution to improve image-based and position-based visual servoings[END_REF][START_REF] Kyrki | New shortest-path approaches to visual servoing[END_REF], then, both the target center-point trajectory in the image and the camera trajectory in the base frame F b will be straight lines. This of course will be true if the system is well modeled. Here, c * t c is the translation vector between F c and F c * ; x o and y o are the image coordinates of target center-point o; θu z is the third component of axis-angle θu which corresponds to the rotation matrix c * R c between F c * and F c . In this paper 2½D VS is used as in [START_REF] Zake | Robust 2½D Visual Servoing of a Cable-Driven Parallel Robot Thanks to Trajectory Tracking[END_REF].

, if s = [ c * t T c x o y o θu z ] T [
Fo Oo Op Fp O b F b Oc Fc base object camera moving-platform A i l i cables B i
Fig. 2. Schematic of a CDPR with a camera mounted on its MP The feature vector s is compared to the desired feature vector s * to compute the error e. To decrease the error e, an exponential decoupled form is selected

ė = -λe (1)
where λ is a positive gain. Camera velocity is then expressed as:

c v c = -λ L -1 s e (2)
where L s is the interaction matrix, defined such that ṡ = L s c v c , and its expression for 2½D VS is given in [START_REF] Chaumette | 2 1/2 D visual servoing: a possible solution to improve image-based and position-based visual servoings[END_REF]. The MP twist, expressed in F p , is defined from the camera velocity c v c as follows:

p v p = A d c v c (3)
where A d is the estimation of the adjoint matrix A d [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]:

A d = p R c [ p t c ] × p R c 0 3 p R c (4)
The Forward Jacobian matrix A of the manipulator relates the cable velocity vector ˙l to the MP twist p v p :

˙l = A p v p (5)
Fig. 2 illustrates the main components of a CDPR and the frames of interest. For the modeling, the cables are supposed to be non-elastic and massless. The pulleys are assumed to be small, thus their radius is neglected. The components of matrix A can then be expressed in the MP frame F p as [START_REF] Pott | Cable-Driven Parallel Robots: Theory and Application[END_REF]:

A =    p u T 1 ( p b 1 × p u 1 ) T . . . . . . p u T m ( p b m × p u m ) T    (6)
where p u i is the unit vector of p # » A i B i , pointing from the exit point A i to the anchor point B i of the ith cable, and p b i is the vector pointing from the origin of MP frame F P to B i .

Vector p u i is known from the MP pose estimation in the base frame F b . As recalled in the introduction and described in [START_REF] Zake | Vision-Based Control and Stability Analysis of a Cable-Driven Parallel Robot[END_REF][START_REF] Zake | Robust 2½D Visual Servoing of a Cable-Driven Parallel Robot Thanks to Trajectory Tracking[END_REF], it is indeed possible to coarsely estimate the MP pose by successive integration of the control output.

The model of the system shown in Fig. 1 is written by injecting Eqs. ( 2), ( 3) and ( 5) into (1):

ė = L s A -1 d A † ˙l (7)
The output of the control scheme, that is, the cable velocity vector ˙l, is expressed by injecting ( 2) and ( 3) into ( 5):

˙l = -λ A A d L -1 s e (8)
l i = 2 m l i = 2 :1 2 m l = 2 m l = 1 :8 8 m A B A B A B A B _ l = :12 m/s for 1 s _ l = :12 m/s for 1 s Fig. 3. Effect of cable slackness on MP displacement
The closed-loop equation of the system is expressed by injecting ( 8) into ( 7):

ė = -λ L s A -1 d A † A A d L -1 s e (9)
According to Lyapunov stability analysis [START_REF] Khalil | Nonlinear systems[END_REF], the stability criterion Π is defined as:

Π = L s A -1 d A † A A d L -1 s ( 10
)
For the system to be stable, Π > 0 is a sufficient condition. Indeed, as long as this condition is fulfilled, it is clear from ( 9) that the error e will converge to 0.

III. TENSION CORRECTION

Due to modeling errors and the estimation of the MP pose, the cables can become slack during task execution. While some limited slackness does not make the system unstable, it acts as an additional perturbation. Indeed, if a cable is slack, applying a velocity to it will not produce the desired displacement of the MP. An example is shown in Fig. 3. Here, on the left a cable is in tension with initial length l i = 2 m and after applying a velocity l = 0.12 m/s for 1 s the final cable length is l f = 1.88 m and the MP has been moved. On the right the initial location of points A and B is the same as before, but the cable is now slack and its initial length is l i = 2.12 m. With the same cable velocity, the final cable length becomes l f = 2 m, but the MP has not moved.

Since with slack cables the actual displacement of the MP differs from the estimated one, the MP pose estimation accumulates additional errors. This in turn leads to computation of cable velocities that may increase cable slackness.

Cable slackness can be detected in several ways. A set of four stereo-cameras is used in [START_REF] Dallej | Vision-based modeling and control of large-dimension cable-driven parallel robots[END_REF] to determine the cable angle, which in combination with MP pose measurement allowed the system to detect cable slackness. Kraus et al. in [START_REF] Kraus | System identification and cable force control for a cable-driven parallel robot with industrial servo drives[END_REF] used tension sensors to directly measure cable tensions. The measurements were then compared to the output of a TDA and a cable length correction was produced. The final control output was a sum of the inverse kinematics output and the correction. From their experiments, all the cables of a fully-constrained CDPR are in tension and the cable tensions have tripled after correction. As a consequence, the stiffness of the manipulator is substantially increased.

Unlike [START_REF] Kraus | System identification and cable force control for a cable-driven parallel robot with industrial servo drives[END_REF] where a position control scheme is developed, our control is in velocity. Thus, it is necessary to relate the cable tension measurements with cable velocities. To do this we propose the control scheme shown in Fig. 4. Here, at every iteration cable tensions are measured and recorded in vector τ m . Each component of this vector is then compared to a threshold tension τ t and the ith component δ li of the cable velocity correction vector δ ˙l is computed as:

δ li = -k c (τ t -τ mi ) if τ mi < τ t , i = 1, . . . , m 0 otherwise ( 11
)
where τ mi is the measured tension of the ith cable, and k c is a positive gain that needs to be tuned. The tension τ t is the value under which the cable is considered to be slack. The resulting behavior depends on the sign of ith cable velocity li computed by the controller. Note that negative velocity leads to the cable becoming shorter, and positive velocity leads to the cable becoming longer. Furthermore, according to (11) δ li is not positive. Thus applying δ li to a negative li results in cable length reducing faster. Applying δ li to a positive li results in cable length increasing slower.

It should be noted that the correction speed is a function of the gain k c . The greater the gain, the faster the correction. However, setting k c too high may perturb the main controller. Thus, k c should be tuned based on the frequency of the TCA loop. Furthermore, k c is used to change the order of magnitude between the tension difference and cable velocities. For example, if (τ t -τ mi ) ≈ 1 N and li ≈ 0.05 = 5 × 10 -2 m/s, then k c will be defined as k c = 10 -2 m/Ns.

The threshold τ t is a tension that is feasible for every cable no matter the MP pose. Such a tension can be found by tracing a workspace, e.g. the Static Feasible Workspace (SFW) with ARACHNIS software [START_REF] Ruiz | ARACHNIS: Analysis of Robots Actuated byCables with Handy and Neat Interface Software[END_REF], and choosing the lower tension bound τ lb . Indeed, τ lb is feasible for every cable and for all MP poses within SFW, thus it is a good pick for τ t .

The proposed tension correction algorithm (TCA) is greatly simplified when compared to [START_REF] Kraus | System identification and cable force control for a cable-driven parallel robot with industrial servo drives[END_REF], while the overall control is more robust thanks to the use of VS. Indeed, instead of using a TDA and a complex computation of the final tension correction, we simply compare the current tensions to a threshold. In this case, the need for the knowledge of the MP pose is avoided since the tension threshold is constant. It is very convenient, because we only have a coarse MP pose estimation. Furthermore, using a coarse MP pose estimation with a TDA could be impossible. More precisely, the tension set provided by the TDA could be unattainable with the current MP pose in case of estimation errors. This would lead to perturbing the VS controller and possibly even making it unstable, thus failing the task. In addition, the simpler the calculations, the faster the controller response to each incoming image and tension measurement. The chosen expression of the CDPR model and thus stability criterion Π does not contain cable tensions. For CDPRs with light cables, if one cable becomes slack, as detected by tension measurements, then it will be considered non-existent. Indeed, a light cable exerts a negligible force on the MP when slack. Furthermore, while a cable is slack, reducing its length does not induce any MP motion. For the stability analysis the expression of the Jacobian matrix A is updated at every iteration accordingly. For example, if at t = 0 s tension τ m1 < τ t , then the corresponding row of matrix A is removed and it becomes a (7 × 6)-matrix. Furthermore, if more than two cables are slack, the Jacobian A will be rank deficient and the robot will be underactuated. Considering the stability criterion Π, it is clear that once rank(A) < 6, then rank(Π) < 6 and thus the stability criterion will not be fulfilled, as shown in Section IV-C.

IV. EXPERIMENTAL VALIDATION

A. Experimental Setup

A six-DoF suspended CDPR prototype with eight cables, named ACROBOT and shown in Fig. 5, is used to validate experimentally the effect of cable slackness on the system behavior. Its frame is a 1.2 m × 1.2 m × 1.2 m cube. The MP size is 0.18 m × 0.17 m × 0.07 m and its mass is 3.5 kg. The cables are Dyneema SK78 2 mm. Tension sensors are mounted on the cables close to their anchor points B i . They are calibrated for a range from -25 N to 75 N, with an accuracy of 0.24 N and repeatability of 0.01 N.

A camera is mounted on the MP facing the ground. The image resolution is set to 640 × 480 pixels. The camera observes one AprilTag [START_REF] Olson | AprilTag: A robust and flexible visual fiducial system[END_REF], a fiducial marker that is used as a target, and algorithms from the ViSP library [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF] recognize and localize it once it becomes visible in the image. From these measurements, the feature vector s is computed and compared to the desired feature vector s * , and the robot is controlled using [START_REF] Chellal | A Kinematic Vision-Based Position Control of a 6-DoF Cable-Driven Parallel Robot[END_REF] to reduce this difference.

The main loop frequency is 25 Hz. The TCA loop can be as fast as 512 Hz, but for the following experiments it was kept at 25 Hz. The tuning of gain λ for 2½D VS on ACROBOT is detailed in [START_REF] Zake | Robust 2½D Visual Servoing of a Cable-Driven Parallel Robot Thanks to Trajectory Tracking[END_REF]. By tracing SFW with ARACHNIS [START_REF] Ruiz | ARACHNIS: Analysis of Robots Actuated byCables with Handy and Neat Interface Software[END_REF], it was found that setting the lower bound tension τ lb = 1 N ensures that all cables are in tension, while keeping the workspace large. Thus, we set τ t = τ lb = 1 N . Finally, TCA gain k c is tuned to be 0.04 m/Ns.

B. Experiments with cable slackness at initial MP pose

For the first set of experiments, a fixed AprilTag, shown in Fig. 5, is observed by the onboard camera. The initial values are the following: Every experiment starts at the initial MP pose b p p0 corresponding to the initial feature vector s. Similarly, the desired MP pose b p * p , corresponding to the desired feature vector s * , is the same for all experiments.

According to [START_REF] Merlet | Simulation of discrete-time controlled cable-driven parallel robots on a trajectory[END_REF] only six cables need to be in tension to achieve any given pose with a suspended CDPR. Given our initial MP pose b p p0 , slack can be introduced on l 1 and l 2 . Once the MP was in its initial pose b p p0 , the winches of cables C 1 and C 2 were turned to increase their length by 4 cm or 8 cm. Six experiments, named E 1 to E 6 , are defined: 6 and7. Please also refer to the attached video to see the recorded experiments. Cable lengths, velocities, and tensions for cables C 1 , C 2 , C 5 and C 6 are shown in Fig. 6. The curves for cables C 3 , C 4 , C 7 and C 8 are not shown. Due to the symmetry in CDPR design and the diagonal MP trajectory, they do not contain any significant information. The added slack can be seen in Figs. 6a and6b, where cable lengths l 1 and l 2 are shown.

• E 1 is 2½D VS without cable slack; • E 2 is 2½D VS
First, the behavior of the CDPR under 2½D VS control with slack and no TCA is analyzed. Cables C 1 and C 2 remain longer in E 3 and E 5 than in E 1 along most of the trajectory (gray and green curves against blue curve, resp., in Figs. 6a and6b). Note that by the end of the trajectory the slack has been transferred to cables C 5 and C 6 , as t, s 0.9

E1: 2½D VS, no slack E2: 2½D VS with TCA, no slack E3: 2½D VS, slack=4 cm E4: 2½D VS with TCA, slack=4 cm E5: 2½D VS, slack=8 cm E6: 2½D VS with TCA, slack=8 cm Fig. 6. Cable lengths l i , cable velocities li and cable tensions τ i for cables C 1 , C 2 , C 5 and C 6 on first, second and third rows, respectively shown in Figs. 6c and6d. Indeed, the six-cable configuration corresponding to the final MP pose is different from the initial one. The transfer of slack can also be seen in the tension measurements. Here, at first τ 1 and τ 2 are below the threshold, but by the end of the trajectory their values increase, while τ 5 and τ 6 decrease below τ t . Note that more than two cables can become slack. For instance, at t = 4 s for both the gray and green curves, the cable tensions τ 2 , τ 5 and τ 6 are below τ t .

Once TCA is added to the controller, the behavior is different. It can be seen that the brown curve of E 4 aligns with the blue curve at approximately t = 0.8 s in Figs. 6a and6b. Thus, at t = 0.8 s cable lengths l 1 and l 2 are equal for blue and brown curves of E 2 and E 4 , resp. This is because up until t = 0.8 s velocities l1 and l2 were increased by the TCA velocity correction, as shown in Figs. 6e and6f. Furthermore, at t = 0.8 s cable tensions τ 1 and τ 2 start to increase. Thus, at this moment cables are no longer slack. In E 6 , shown in cyan, the correction of 8 cm slack takes about twice the time. Velocities l1 and l2 are increased until approximately t = 1.5 s, when the tensions τ 1 and τ 2 start to increase and cable lengths l 1 and l 2 align with the blue curve. Finally, the velocity and length curves for cables C 5 and C 6 are the same during experiments E 2 , E 4 and E 6 , meaning that there was no need for slack correction after the change of six-cable configuration.

The AprilTag center-point trajectory in the image is visible in Fig. 7. The larger the initial cable slack, the larger the deviation from the straight-line trajectory. As the MP rapidly falls to a new cable configuration, it causes a sideways drift in the image trajectory. On the contrary, when TCA is used, the produced trajectories are very close to the desired straight lines. In conclusion, TCA improves the VS controller by reducing cable slack. Eventually, cable slackness affects the accuracy of MP pose estimation, because the desired motion is not produced. In Fig. 8 the distance between the measured final MP pose and its estimation is shown for each experiment. Here, the rotational deviation is defined as the angle θ of the axis-angle θu representation of the rotation matrix p R p between the estimated and the actual MP pose. In E 1 , E 2 , E 4 and E 6 the estimation errors amount to 0.03 m and 2 • . Thus, using TCA improves the MP pose estimation. It is significantly better in case of initial slack when compared to the classic controller. In E 5 the errors reach 9 cm and 5 • . Given the short trajectory time of only 8 s, it can be concluded that large cable slack leads to fast deviation of the MP pose estimation. Moreover, the larger the initial cable slack, the larger the error in MP pose estimation.

C. Stability Analysis

This section deals with the stability analysis for the given experiments. The stability criterion Π is computed off-line as defined in [START_REF] Ramadour | Grasping Objects With a Cable-Driven Parallel Robot Designed for Transfer Operation by Visual Servoing[END_REF] MP via a static camera shown in Fig. 5. Furthermore, cable slackness is also taken into account as described at the end of Section III.

Fig. 9 contains the results of stability analysis over time. Here, the boolean Stable is true if Π > 0 and false otherwise. As expected, when no cable slack is added to the system, it remains stable throughout the trajectory no matter the controller, as shown in Fig. 9a. At one iteration one tension measurement goes below the threshold, however that does not make the system unstable, because there are still seven cables in tension. For E 3 , shown in Fig. 9b in gray, at t = 2.5 s three cables become slack. Furthermore, at t = 3.2 s even four cables become slack. This leads to Π > 0 no longer being true. Thus, when the MP is underactuated, it is no longer in static equilibrium and the stability of the system is not ensured. Hence, the static equilibrium of the MP and the stability of the system are linked. However, if TCA is used, it rapidly corrects the slack and the system remains stable all along the trajectory (see brown curve in Fig. 9b).

The behavior is analogous in E 5 and E 6 , shown in Fig. 9c. Due to larger initial slack in E 5 , the time period of underactuation is increased when compared to E 3 . Similarly to E 4 , in E 6 TCA successfully reduces slack, thus avoiding underactuation of the MP and loss of stability.

D. Tracking a mobile robot

The attached video contains a second experiment dealing with the MP following a mobile robot moving on a track over a prolonged time period. The velocity of the mobile robot is unknown, the VS controller is simply tasked to keep its AprilTag in the desired pose. To aid with MP pose estimation, a known tag is put on the track and the estimation is corrected every time it becomes visible. Despite this, the classic VS fails due to cable slack in lap 29 at t < 10 min. On the contrary, VS with TCA was perfectly functional after 700 laps, at t > 180 min.

V. CONCLUSIONS

This paper proposed the use of a tension feedback loop to deal with cable slackness. The algorithm is simple and improves the visual servoing controller. Indeed, the tension correction only occurs when cable slackness is detected and it is stopped as soon as the tension threshold is reached. The proposed algorithm is efficient in reducing cable slack and its rapidity depends on the tuning of the gain and the amount of slack on the cable.

Cable slack affects not only the CDPR responsiveness to control output, but also its stability. Indeed, when initial slack is large, the MP can become underactuated along the trajectory and that leads to the stability criterion being no longer ensured. In addition, the larger the cable slack, the larger the MP pose estimation error, when TCA is not used.

When cable slack is transferred between cables, the MP can sway uncontrollably, leading to sharp changes of the target trajectory in the image. As the cable slack increases, it becomes more likely that the target will leave the field of view, resulting in a task failure.

Future work includes the validation of the proposed controller on a fully-constrained CDPR and on a large-scale CDPR with elastic cables.

Fig. 1 .

 1 Fig. 1. Control scheme for visual servoing of a CDPR

Fig. 4 .

 4 Fig. 4. Visual servoing control of a CDPR with a TCA

Fig. 5 .

 5 Fig. 5. ACROBOT: a CDPR prototype located at IRT Jules Verne, Nantes

•

  b p p0 = 0.202 m; 0.118 m; 0.268 m; -18 • ; 10 • ; 4 • • c p o0 = 0.1 m; -0.05 m; 0.26 m; 165 • ; 10 • ; 179 • • o 0 = 0.38 m; -0.19 m and desired values are selected to be: • b p * p = -0.11 m; -0.20 m; 0.366 m; 13 • ; -20 • ; 33 • • c p * o = -0.14 m; 0.115 m; 0.35 m; 178 • ; -20 • ; 147 • • o * = -0.39 m; 0.33 m where b p p is the MP pose in F b ; c p o is the AprilTag pose in F c ; and o stands for the AprilTag center-point coordinates in the image. The angles in pose parameters are Euler angles. Note that c p oi and o i are measured at every iteration i to be used in the control scheme. Meanwhile, b p p0 and b p * p have been obtained from an external camera observing another AprilTag mounted on the MP, shown in Fig. 5, and serve as ground truth. It is different from the b p p used in the control.

Fig. 8 .

 8 Fig. 8. Translational and rotational deviation of final MP pose estimation

Fig. 9 .

 9 Fig. 9. The relation between the amount of cables in tension and system stability for: (a) no slack; (b) 4 cm slack; (c) 8 cm slack

TABLE I NOTATION

 I USED THROUGHOUT THE PAPER

  with TCA without cable slack; • E 3 is 2½D VS with 4 cm cable slack; • E 4 is 2½D VS with TCA with 4 cm cable slack; • E 5 is 2½D VS with 8 cm cable slack; • E 6 is 2½D VS with TCA with 8 cm cable slack. The experimental results are shown in Figs.

  using the following recorded variables:

	Initial position			E5: 2½D VS, slack=8 cm
	Desired position		E2: 2½D VS with TCA, no slack
	E1: 2½D VS, no slack		E4: 2½D VS with TCA, slack=4 cm
	E3: 2½D VS, slack=4 cm	E6: 2½D VS with TCA, slack=8 cm
		0
		100
	v, pixels	200 300
		400
			0	100	200	300	400	500	600
					u, pixels
	Fig. 7. AprilTag center-point trajectory in the image
	0 2 4 6 10 8	Translational deviation, cm Rotational deviation, degrees
			E1	E2	E3		E4	E5	E6

This work is supported by IRT Jules Verne (French Institute in Research and Technology in Advanced Manufacturing Technologies for Composite, Metallic and Hybrid Structures) in the framework of the PERFORM project.


Utkarsh A Mishra
email: umishra@me.iitr.ac.in

Marceau Métillon
email: marceau.metillon@ls2n.fr

Stéphane Caro
email: stephane.caro@ls2n.fr

Ieee Member

Kinematic Stability based AFG-RRT* Path Planning for Cable-Driven Parallel Robots †

Motion planning for Cable-Driven Parallel Robots (CDPRs) is a challenging task due to various restrictions on cable tensions, collisions and obstacle avoidance. The presented work aims at proposing an optimal path planning strategy in order to both maximize the wrench capability and the dexterity of the robot in a cluttered environment. First, an asymptoticallyoptimal path finding method based on a variant of rapidly exploring random trees (RRT) is implemented along with the GilbertJohnsonKeerthi (GJK) algorithm to account for the collision detections. Then, a goal biased Artificial Field Guide (AFG) is employed to reduce convergence time and ensure directional exploration. Finally, a post-processing algorithm is added to get a short and smooth resultant path by fitting appropriate splines. The proposed path planning strategy is analyzed and demonstrated on a simulated and experimental setup of a six-DOF spatial CDPR.

I. INTRODUCTION

Cable-Driven Parallel Robots (CDPRs) belong to the class of parallel robots in which the end-effector (EE), a movingplatform, is connected to a rigid frame by means of cables instead of articulated limbs. CDPRs can have a low inertia, a higher payload to weight ratio and a large workspace which enables a variety of usage of such robots in large scale 3D printing [START_REF] Barnett | Large-scale 3D printing with a cable-suspended robot[END_REF], large scale telescopes [START_REF] Zi | Dynamic modeling and active control of a cable-suspended parallel robot[END_REF] and in rehabilitation mechanisms [START_REF] Homma | A wire-driven leg rehabilitation system: Development of a 4-DOF experimental system[END_REF]. However, there are some challenges with such systems along with the advantages accompanying it. The cables used in such robots can only apply unilateral forces, i.e. they can only pull the moving platform and not push it [START_REF] Gouttefarde | Intervalanalysis-based determination of the wrench-feasible workspace of parallel cable-driven robots[END_REF].

Given that, the designing of CDPRs based on manipulability and structural stability to external wrenches plays a vital role in governing the workspace and reliability under varying conditions. Gosselin and Angeles first used the condition number of the kinematic Jacobian matrix in 1990 [START_REF] Gosselin | Singularity Analysis of Closed-Loop Kinematic Chains[END_REF] [START_REF] Gosselin | Global performance index for the kinematic optimization of robotic manipulators[END_REF] which was extended by Cui and Hao [START_REF] Cui | Kinematic performance analysis on a new spatial rotation 3-DOFs parallel robot mechanism[END_REF] as a measure of the mechanism dexterity. The problem of non homogeneous jacobian was addressed by Guay et al. [START_REF] Guay | Measuring how well a structure supports varying external wrenches[END_REF] which further explored an index used for supporting varying force and motion capabilities of the actuator, called the capacity margin. A similar index [START_REF] Lessanibahri | Twist feasibility analysis of cable-driven parallel robots[END_REF] was also established for twist feasibility based on cable velocities. Various works [START_REF] Rosati | On the design of adaptive cable-driven systems[END_REF], [START_REF] He | Kinematic performance evaluation of 2-degree-of-freedom parallel mechanism applied in multilayer garage[END_REF] have used related indices as a measure of the effectiveness of the Based on a particular configuration, motion planning for CDPRs has been explored by means of various strategies. Cluttered environment with obstacles introduces many challenges related to obstacle avoidance, cable collision and management of cable tensions. A dynamic point-to-point trajectory planning for underactuated CDPRs was proposed in [START_REF] Id | Rest-to-Rest Trajectory Planning for Underactuated Cable-Driven Parallel Robots[END_REF]. Previously, various strategies to plan point-to-point trajectories using polynomial and trigonometric functions were introduced in [START_REF] Gosselin | Experimental determination of the accuracy of a three-Dof cable-suspended parallel robot performing dynamic trajectories[END_REF], [START_REF] Gosselin | Dynamic point-to-point trajectory planning of a two-DOF cable-suspended parallel robot[END_REF] for navigating through waypoints. Point-to-point motion planning is required for pickand-place like operations. However, these works only validated the wrench feasibility conditions without obstacle avoidance. Motion planning of reconfigurable Mobile CD-PRs with Direct Transcription Method while considering the static equilibrium and the tipping conditions of the overall system along the searched path was further explored [START_REF] Rasheed | Optimization based Trajectory Planning of Mobile Cable-Driven Parallel Robots[END_REF] along with obstacle avoidance.

In context of path planning, sampling based path planning algorithms like rapidly-exploring random trees (RRT) perform effective planning in cluttered environments. RRT*, a variant of general RRT, will guarantee asymptotic optimality [START_REF] Karaman | Sampling-based algorithms for optimal motion planning[END_REF] in finding path to a goal, if it exists. A modified RRT* along with a fast GJK obstacle avoidance algorithm [START_REF] Gilbert | A Fast Procedure for Computing the Distance Between Complex Objects in Three-Dimensional Space[END_REF] [18] [START_REF] Karaman | Sampling-based algorithms for optimal motion planning[END_REF] was proposed in [START_REF] Bak | Collision-free path planning of cable-driven parallel robots in cluttered environments[END_REF] considering cable/cable and cable/moving-platform collisions as well as the wrench capability of the system. A faster convergence criterion based on the condition number of the kinematic Jacobian matrix with the help of an artificial potential field was introduced in [START_REF] Zhang | Optimal RRT* Planning and Synchronous Control of Cable-Driven Parallel Robots[END_REF]. Recently, RRT* has also been used to compute efficient trajectories based on horizontal linear periodic motion energies of the state [START_REF] Xiang | Dynamic Point-To-Point Trajectory Planning for Three Degrees-of-Freedom Cable-Suspended Parallel Robots Using Rapidly Exploring Random Tree Search[END_REF].

In this paper, RRT* based offline path planning has been explored for pick-and-place operations (point-to-point motions) in a cluttered environment. A pose based capacity margin as a measure of robustness on sustaining external wrenches is used as the cost function in addition to the dexterity. Finally, the AFG-RRT* algorithm is formulated and solved under the guidance of a goal-biased Artificial Field Guide (AFG) and a GJK based detection of cable/cable, cable/obstacles and moving-platform/obstacles collisions. The overall objective is to plan the path to the goal, if it exists, maximizing the structure stability to external wrenches and with proper manipulability. The motion was simulated and demonstrated on a 6-Degrees Of Freedom (DOF) spatial CDPR shown in Fig. 1.

Section II describes the robot model and presents the problem statement. Section III formulates the cost indices. Section IV deals with the implementation of the AFG-RRT* algorithm along with the functions of transition between poses and collision detection. The simulated and experimental results are discussed in Section V. Finally, conclusions are drawn and future work is presented in Section VI.

b a i + l b i u i -b R p p b i -b p = 0 3 ( 1
)
u i is the unit vector along the ith cable pointing from A i to B i . The subscript b ( p , resp.) means that the corresponding vector is expressed in frame F b (F p , resp.).

The upper and lower bounds on cable tensions are denoted as t and t, respectively (t = 86 N, t = 1 N). The MP mass is named m E .

This paper aims to define a methodology to generate optimal paths in a cluttered environment from a given starting pose to a given ending pose while maximizing the overall robot performances. The maximum and minimum tension impose significant restrictions to the motion of the moving platform. The only external wrench action on the moving platform is due to its own weight. The practical setup shown in Fig. 1 of size 4.24 m x 3.67 m x 2.76 m (l x b x h) has a 5.6 Kg moving platform of size 0.28 m x 0.28 m x 0.20 m.

III. PERFORMANCE INDEX A. Capacity Margin Formulation

Let an external wrench, w e , is applied on the moving platform given by, w e = [f T p τ τ τ T p ] T for an external force f p and moment τ τ τ p . For the moving platform to remain in equilibrium, there must exist a m-dimensional vector of cable tensions, t = [t 1 ,t 2 , . . . ,t m ], in order to satisfy the Newton-Euler equation given by,

-Wt + w e = 0 6 (2)
where, W is the Wrench matrix of the mechanism at this particular pose of the moving platform, given by

W = b u 1 . . . b u m (1/r) b R p p b 1 × b u 1 . . . (1/r) b R p p b m × b u m (3)
and

w e = [f T p (1/r)τ τ τ T p ]
T reformulated as dimensionally homogenized matrices with the help of a characteristic length r. This length is determined by imposing isotropy conditions on the wrench matrix. For planar manipulators [START_REF] Rosati | On the design of adaptive cable-driven systems[END_REF], it is given by

r 2 = (2/m) ∑ m i=1 F b i 2 2
, whereas, for spatial manipulators, it is obtained for minimum Euclidean distance between square matrix WW T and diagonal matrix of squared eigenvalues (σ 1 , . . . , σ n ) of WW T .

Given the above formulation, the equilibrium in eq. ( 2) must hold for all the possible external wrenches (w e ) which can be specified as path planning requirements. Let us denote this set of external wrenches as W e in R 6 , given by W e = {w e ∈ R 6 } represented as a polytope in the wrench space of the robot. Similarly, due to bounded cable tensions, let the allowable tension set, T be composed of all the feasible cable tension vectors such that

T = {t ∈ R m : t ≤ t ≤ t} (4)
represents another polytope (a hypercube) in the tension space of the robot. The two sets W e and T are defined in two different spaces, namely the wrench space and the tension space, which are connected by eq. ( 2). The feasible tensions t ∈ T are linear mapped onto the wrench space using the transformation wrench matrix, W, to get the set of available feasible wrenches, w t . Hence, after transformation, T becomes a special type of polytope in wrench space called the zonotope, given by W t = {w t :

w t = Wt ∀ t ∈ T } (5)
and referred to as the Available Wrench Set (AWS). This can be visualized from Fig. 3 as a low dimensional representation of tension space (3 dim.) and wrench space (2 dim.).

Fig. 3: Cable tension space, Available wrench set, External wrench set and capacity margin index (γ C )

For the system to be able to reach an equilibrium, there must be available cable tensions to balance the external wrench, i.e.,

W e ⊆ W t (6)
This can be analysed graphically as the degree of inclusion of W e in W t i.e. the distance between two polytopes which is represented as a signed distance, s, that is positive when W e is completely within W t , zero when they touch each other and negative when W e is partly outside W t . Mathematically, the Hyperplane Shifting Method (HSM) [START_REF] Guay | Measuring how well a structure supports varying external wrenches[END_REF] is implemented to convert the zonotope, W t , to a set of hyperplanes satisfying Ax = b, where x ∈ R 6 is a vector in wrench space. Finally, between the ith hyperplane, a i x = b i , and j th task wrench, w T e j , the signed distance, s i j , is given as

s i j = (b i -a i w T e j ) a i (7)
and capacity margin (γ C ) for such a pose is defined as the minimum degree of such inclusion, given by

γ C = min (s i j ) (8)

B. Dexterity Margin Formulation

Let the linear and angular velocities of the moving platform in Cartesian space be denoted as b ν and b ω ω ω respectively. The Kinematic Jacobian matrix, J, defined as

b ν b ω ω ω = J θ θ θ (9)
linearly maps the actuated cable velocities θ θ θ to the twist for a particular pose and by calculation, J = W T . While analyzing transfer of error for the transformation from cable space to the cartesian space in eq. ( 9), effective manipulability can be well ensured by reducing the magnitude of the transfer of error and can be related to the dexterity of the robot [START_REF] Merlet | Jacobian, manipulability, condition number, and accuracy of parallel robots[END_REF], which is the second performance index considered in the paper. The conditioning number of the kinematic Jacobian matrix is defined as,

κ = λ max (J T J) λ min (J T J) (10)
where, λ max , λ min denote the maximum and minimum eigenvalues of matrix J T J, σ max , σ min denote the maximum and minimum singular value of the matrix. The Jacobian matrix conditioning number, κ is a measurement for the dexterity and it's a number no less than 1, i.e. 1 ≤ κ ≤ ∞. Hence, the dexterity margin, γ D , is formulated as the reciprocal of this conditioning number i.e.

γ D = 1/κ ( 11
)
and thus 0 ≤ γ D ≤ 1. When γ D = 1, the mechanism is at its optimal transmission capacity. Whereas, when γ D = 0 the mechanism loses its control and attains a singular configuration. Thus γ D is used to measure a robots dexterity and controlling precision.

IV. MODIFIED AFG-RRT* FORMULATION A. RRT* Algorithm

The implemented RRT* algorithm, as in Algorithm 1, includes the cost of traveling to each vertex relative to its parent vertex based on the cost function which is modelled as a combination of the capacity, γ C , and dexterity, γ D , margins defined in the previous section. A simple path

Algorithm 1: Modified RRT* Result: Path from Start to Goal Node while True do for V in G do if Distance(V, Goal) < (ε, δ ) then return G; end end X rand = RandomNode(AFG, Obstacle); X nearest = Nearest(G, X rand ); X new = Transition(X rand , X nearest ); if PATH(X rand , X nearest ) ∈ Obstacle then Continue; else X new .Cost = Cost(X new , X nearest ); X best , X neighbors = Neighbors(G, X new , D N ); Link = Chain(X new , X best ); for X ∈ X neighbors do if X new .Cost + Cost(X new , X) < X.Cost then X.Cost = X new .Cost+Cost(X new , X); Parent(X) = X new ; G += {X new , X}; end end G += Link; end end
cost is the cumulative increment in capacity margin and dexterity achieved during the transition. Points are generated over the workspace and connected to the closest available node based on the cost metric. For a graph, G(V, E), of edges and vertices initialized with a start node, the closest node (X new ) is obtained. Now, if a node with a better cost (X best ) is obtained within a neighborhood of existing vertices (X neighbours ) marked by a fixed distance radius, D N , it will replace the proximal node using the distance metric. The algorithm ends when a node is generated within the goal region bounded by ε Cartesian distance and δ angular distance. The overall transition between intermediate nodes in governed by a pose transition rule and checked for obstacle avoidance.

B. Collision Detection

In order to account for the large number of collision checks in sampling based algorithms like RRT*, the GJK algorithm is integrated for its faster and more accurate collision detection behavior. It calculates the distance between objects, with reference to distance from origin. To deal with the limitation of the algorithm to convex polytopes, the moving-platform, obstacles and the cables are approximated as convex objects [START_REF] Bak | Collision-free path planning of cable-driven parallel robots in cluttered environments[END_REF]. All possible interference between the moving-platform, cables and obstacles (shown in Fig. 4) are considered to determine the paths free of collision between two nodes.

p N = p 1 + λ (p 2 -p 1 ) (12)
Q N = sin[(1 -λ )Ω] sin Ω Q 1 + sin[λ Ω] sin Ω Q 2 ( 13
)
where, cos Ω = Q 1 .Q 2 , i.e. the 4-dimensionsal dot product between the two quaternions of the form Q = [q 0 , q1 (3×1) ].

D. Artificial Field Guide (AFG)

Simple uniform random number generators reaches asymptotic optimality with RRT*. Such numbers encourage exploration at the expense of increasing computation time [START_REF] Zhang | Optimal RRT* Planning and Synchronous Control of Cable-Driven Parallel Robots[END_REF]. Thus, a goal biased artificial field guided formulation of generating random numbers is used to attract the search towards the goal and, ensure faster convergence. This way the searching and exploration is much more directional towards the desired goals. Nodes, N(p, Q) are generated based on the AFG, 15)

N(p) = C(p) + (Ds(ψ)c(φ ), Ds(ψ)s(φ ), Dc(ψ)) (14)
P(D | D ∈ [D min , D max ]) ∝ 1 N(p) -C(p) 2 (
N(Q) = sin[(1 -λ r )Ω r ] sin Ω r C(Q) (16)
where s(ψ) = sin ψ and c(ψ) = cos ψ. The probability of choosing D in eq. ( 15) governs the attraction phenomenon where D min is a decision variable and D max is the maximum radius coverage of the AFG (in Fig. 5) extended to D max = ∆X Workspace max . N(p, Q) must be within the wrench feasible workspace (i.e. γ C | @N(p,Q) ≥ 0) and outside any obstacle. Further chosen conditions on decision variables are ψ, φ , Ω r ∈ Random(0, 2π) and λ r ∈ Random(0, 1). The probability of C being the goal node or any random node is defined using a probability P g , i.e., P(C = X goal ) = P g and P(C = X rand ) = 1 -P g .

The AFG, as visualized from Fig. 5, is modelled as a potential field where the random node generation is based on an attraction phenomenon by the goal node X goal . The field also preserves randomness by generating random nodes and the degree of randomness is controlled by adjusting P g . The generated node will be neglected if it lies on the obstacles.

E. Post Processing

The intermediate nodes between the start and the goal node, as obtained by the modified RRT* undergo a postprocessing strategy while retracing the nodes from the goal to the start. Each subsequent node will be connected to the most furthest node, if the path connecting them lies out of the obstacle region. The final output is provided as the post processed output and a spline is fitted based on the final nodes to get a smoother path.

V. RESULTS AND DISCUSSIONS

A. Simulation Results

In this section, the above formulated AFG-RRT* algorithm is simulated in a cluttered environment to obtain near optimal path for point to point motions. The complete strategy is simulated on a 6-DOF spatial CDPR shown in Fig. 1. The rotation of the MP about x and y axes is restricted in order to simulate pick-and-place operations. Therefore, a node for RRT* is defined by a 5-dimensional vector composed of [p, Q] where p = [x, y, z] and Q = [q 0 , q 1 ] which forms a unit quaternion q = q 0 + 0 î + 0 ĵ + q 1 k representing MP position and orientation about z-axis (also denoted as β in Euler angle notation in paper).

The wrench matrix was calculated for each pose without considering orientation i.e. β = 0 rad and the tension space was mapped into the wrench space to get W t . W e is considered to be composed of the gravity wrench of the MP, i.e., w E = [0, 0, m E g, 0] where, g = 9.81 m/s 2 is the gravitational acceleration.

Fig. 6a depicts the isocontours of the capacity margin γ C for a constant and a null orientation (β = 0) of the MP throughout the xz plane and for the mid-plane along ydirection. This margin was normalized by considering its

γ = γ C + γ D ∈ (0, 2) ( 17
)
This metric is used to conduct the AFG-RRT* simulation while analyzing separately for the individual metrices (γ C and γ D ). All the simulations are performed using ©MATLAB with CPU computations on an Intel ®i7-7500U CPU@2.70GHz. The AFG-RRT* was initialized using the parameter values given in Tab. I and the task was formulated The obtained path, orientation about z-axis and the rewinded path according to the post-processing algorithm is given in Fig. 8a, Fig. 8b and Fig. 9 respectively. The overall path is subjected to collision avoidance between the cables, MP and the obstacles. The algorithm was simulated for a number of times and the performance index achieved was recorded along with the number of iterations of AFG-RRT* and the analysis results can be visualized in Tab. II. 215 iterations was required to achieve a mean capacity margin equal to 0.55 and mean dexterity margin of 0.69. The average computation time required for simulating the complete AFG-RRT* simulation was 13 minutes.

B. Experimental Demonstration

Cubical boxes of size 44.5 cm x 44.5 cm x 44.5 cm are arranged in order to imitate the simulated environment with cuboidal obstacles and the CDPR. The environment ) was completed in 60 secs and the mean capacity and dexterity margin along the path were 0.5396 (Simulated: 0.5232), and 0.6355 (Simulated: 0.6345) respectively. The closest proximity between cables and obstacles was obtained as 0.8 cm (Simulated: 1.2 cm), 5.7 cm (Simulated: 8 cm) for two particular instants and that between moving-platform and obstacles was found to be approximately 8 cm (Simulated: 10 cm). The differences being primarily due to the pulleys that have not been considered here in the robot model.

VI. CONCLUSIONS

A structural robustness based path planning was performed based on a goal-biased AFG coupled random sampling RRT* algorithm. The algorithm ensures directional exploration, maximizes kinematic stability, maintains efficient manipulability and achieves fast convergence. The algorithm was simulated and demonstrated for a 6-DOF CDPR setup with 8 cables. An analysis of various performance metrics showed the effectiveness of the proposed methodology in planning a path for a suspended CDPR in a cluttered environment.

Fig. 1 :

 1 Fig. 1: 6-DOF spatial CDPR in a cluttered environment

  II. ROBOT MODEL Let us consider an n-DOF CDPR with m cables. Its ith closed loop is represented in Fig. 2. The frame F b of origin O is attached to the base. The frame F p of origin P is attached to the Moving-Platform (MP). a i denotes the Cartesian coordinates vector of exit point A i expressed in F b . b i denotes the Cartesian coordinates vector of anchor point B i expressed in F p . The MP pose is defined by the vector p pointing from O to P and the unit quaternion b Q p ∈ SU(2) (and equivalent rotation matrix b R p ∈ SO(3)) from F b to F p .

Fig. 2 :

 2 Fig. 2: ith closed-loop of the CDPR

Fig. 4 :

 4 Fig. 4: Convex shape representation of the Moving Platform (MP), Cables (A i B i ) and Obstacles (O i ) for satisfying the prerequisites of GJK Algorithm

Fig. 5 :

 5 Fig. 5: AFG within the complete workspace of the CDPR

Fig. 6 :

 6 Fig. 6: Normalized Cost functions in 4-DOF Spatial CDPR AFG-RRT* Simulations

Fig. 7 :

 7 Fig. 7: Start and end MP poses of the CDPR working in a cluttered environment and boundary of its wrench feasible workspace corresponding to γ C = 0

Fig. 8 :

 8 Fig. 8: Simulated Environment Results

Fig. 9 :

 9 Fig. 9: Smooth path after post-processing

Fig. 10 :

 10 Fig. 10: CDPR prototype and cluttered environment made up of 18 boxes

TABLE II :

 II Performance metrics obtained from 100 simulations of AFG-RRT* algorithm

	Metric	Mean Value	Max Value Standard Deviation
	γ C	0.5542	0.7122	0.1104
	γ D	0.6934	0.8237	0.0978
	γ	1.0941	1.4277	0.1507
	No. of Iterations	215	317	22
	Computing time	13.78 min	15.38 min	1.70 min

https://www.umishra.live/redirects/icra2021.html

This work was supported by the ANR CRAFT project, grant ANR-18-CE10-0004, https://anr.fr/Project-ANR-18-CE10-0004


Yousr Skhiri

Philipp Gruner

Benoît Semin

Quentin Brosseau

Deniz Pekin

Linas Mazutis

Victoire Goust

Felix Kleinschmidt

Abdeslam El Harrak

J Brian Hutchison

Estelle Mayot

Jean-Franc ¸ois Bartolo

Andrew D Griffiths

Valérie Taly

Jean-Christophe Baret

Dynamics of molecular transport by surfactants in emulsions

We consider the dynamics of equilibration of the chemical potential of a fluorophore in a monodisperse emulsion containing droplets with two initially different concentrations of the fluorophore. Although the exchange mechanism involves a single time-scale at the droplet (microscopic) level, the organisation of the droplets determines the exchange dynamics at the population (macroscopic) level. The micelle concentration in the continuous phase and the chemistry of the fluorophore control the microscopic exchange rate while the disorder of the initial condition determines the power-law of the long-time scale, recovered in a minimal analytical model. We also show here that an additive in the droplet such as Bovine Serum Albumine (BSA) acts on the microscopic exchange rate and slows down the exchange process, by increasing the solubility of the fluorophore in the dispersed phase rather than creating a viscoelastic layer at the droplet interface.

Introduction

Droplet-based microfluidics has proven a very effective tool for the miniaturisation and automation of biological assays, for single cell analysis, DNA screening or drug screening [1][2][3] . In order to reliably function, the droplets have to fulfill three major conditions: they must be stable against coalescence, biocompatible, and their components must remain encapsulated over time, three properties controlled by the surfactant molecules (see 4,5 and reference therein). Understanding the transport of molecules between dispersed microcontainers is not only beneficial for the optimization of microreactors for biochemical applications 1 , but is also of relevance to other encapsulation systems which -for example -rely on controlled release of the encapsulated material 6 such as drug delivery systems 7 , for applications such as emulsion polymerisation, or for cell physics 8 . The stability of droplets is an old problem linked to the ageing of emulsions 9,10 . In the case of a cream, a mayonnaise or a pharmaceutical formulation, all droplets of the sample contain the same material and ageing is linked to the transport of the dispersed phase through the continuous phase or through the bilayers that can form between droplets 11 . With the recent applications of droplet emulsions for biotechnological applications, the paradigm for emulsion ageing has slightly shifted. In the case of storage of a chemical compound library or a library of single genes in droplets, all droplets contain different molecules, possibly at different concentrations [12][13][14][15][16] . It becomes then important to understand and control the ageing of such emulsions by the exchange of solute.

The fine control of droplet actuation and of droplet size distribution, and the accessibility of short timescales (∼ ms) in microfluidic emulsification -difficult to achieve in bulk-emulsification -makes it appealing as a new tool to quantitatively study the physics and physical chemistry of interfaces and emulsions [17][18][19][20][21][22][23] . Studies on solute exchange have already been performed in bulk emulsions 9,[24][25][26] , but these studies do not capture the microscopic details of the exchange at the single droplet or interface level. Such insights are now accessible by the microfluidic control of droplets 11,21,27 , but the understanding of collective dynamics in population of droplets is still lacking 28 .

In the present manuscript, we address the dynamics of molecular exchange between droplets in a mixed emulsion. We measure experimentally the relaxation of concentration differences in an emulsion (produced in microfluidics) initially containing droplets with two different concentrations of a fluorescent dye and make the link between the microscopic exchange between two adjacent droplets and the macroscopic kinetics (at the scale of the emulsion). We study the case of water in fluorinated oil inverse emulsion as they represent the most promising systems for biotechnology applications 5 . Fluorinated oils are used to circumvent the exchange of molecules between droplets since organic compounds are highly insoluble in these oils 29 : the chemical stability of encapsulation is improved compared to mineral or organic oils.

Materials and methods

Chemicals

Resorufin sodium salt, fluorescein sodium salt, rhodamine 101 inner salt, sulforhodamine 101 were purchased from Sigma-Aldrich. Solution were prepared by dilution in Phosphate Buffer Saline

Surfactant

We used droplets of water in a fluorinated oil (HFE-7500, 3M) stabilised against coalescence by a block-copolymer surfactant, see Figure 1(a). The surfactant was prepared from the commercially available carboxylic acid Krytox 157-FSH (Dupont) and Jeffamine polyetheramines (ED 600, Huntsmann). We used surfactant solution prepared from two different batches and refered in the following as KryJeffa (1) and KryJeffa (2) . The synthesis scheme was based on the synthesis described in Holtze et al. 30 .

Surfactant characterizations

To characterize our system, the surface tension at different surfactant concentrations were measured for the KryJeffa -HFE mixtures using pendant droplet method (DataPhysics OCA). Equilibria are reached after a time of the order of several minutes to several hours, depending on surfactant concentration for droplets of volumes of ∼ 1 to 20 µL. Using Dynamics Light Scattering (Malvern Zetasizer), we measured the scattered intensity of the surfactant solution as a function of concentration. We used the automatized optimised mode for signal acquisition and rescaled the data by the corresponding gain provided by the manufacturer. From the scattered intensity, we deduced the Critical Micellar Concentration (CMC) of the surfactant. The size of the scattering objects was obtained by the time-correlation of the scattered signal. The two methods lead to a value of CMC of the surfactant in HFE-7500 in the range 0.02-0.04% (Figure 1(b,c)). The DLS measurements also showed an increase of the scattering signal corresponding to an increase of the number of objects with the surfactant concentration (quasi-linear relationship) with a weak change in the size of the reverse micelles at higher concentrations (from a radius of 80 nm to 120 nm over 3 decades in concentration).

Microfabrication

Microchannels of depth 60 µm were produced with standard soft-lithography techniques of replica molding of a SU-8 master in Polydimethylsiloxane 31 . The PDMS stamp was then bound to glass after oxygen plasma activation and hydrophobized using a commercial coating agent (Aquapel, PPG Industries). A collection vial was prepared using a 1 mL plastic syringe (BBraun) cut at both ends with a scalpel. Two cylindrical PDMS blocks cut from a crosslinked PDMS slab were prepared: in one PDMS cylinder, one hole was punched with a Harris Unicore biopsy punch (0.75 mm diameter) to allow a connection with a teflon tubing. Two holes were punched in the second cylinder. Both cylinders were then used to close both ends of the vial.

Chip connection

All flow rates were controlled by Nemesys syringe pumps (Cetoni GmbH). The bottom of the collection vial was connected to the outlet of the droplet production device at one hole and to a syringe at the second one. The syringe was actuated by the syringe pumps to either push or withdraw liquid (Control Pump). The top of the collection vial was connected back to the chip. In a preliminary phase the collection vial was filled with oil-surfactant solution (Figure 2 (a)).

Droplet production, collection and reinjection

Droplet-based microfluidic systems were used to produce monodisperse emulsion. Droplets of resorufin in PBS solutions (0 µM and 10 µM) were produced in parallel at two different nozzles by flow focusing in a fluorinated oil (HFE-7500, 3M). The droplets contained a fluorescein concentration used to encode the droplet type (150 nM and 50 nM). The measurement of the fluorescein fluorescence provided a way to determine if droplets initially contained 0 µM or 10 µM resorufin, even after the resorufin has equilibrated and the concentration in all droplets is identical.

The oil flow rate was 100 µL/min (except when BSA is used, the oil flow rate being increased to 200 µL/min), and each acqueous flow rate 50 µL/min corresponding to droplet volume of 0.3 nL.

Droplets were stabilised against coalescence by the surfactant at concentrations ranging from 0.1% (limit of emulsion stability) to 5% (w/w). No coalescence was observed on chip. As demonstrated in a previous study 32 , preventing droplets from touching each other after generation increases emulsion stability. Therefore we designed the channels with an incubation line of about 1 centimeter to increase droplet stability. During the stabilization time of the whole system, the Control Pump CP (Figure 2(b)) was set up to 'infuse' and the droplet flow towards the outlet. When the system was stabilised, CP was set to 'withdraw' with a flow rate tuned to stop the flow at the outlet of the collection vial. The emulsion was collected in the vial. After a few minutes (typically 5 minutes), a volume of about 500 microlitres of emulsion was stored in the vial (∼ 10 6 droplets).

The production of the droplets was stopped, CP was set to 'infuse' (typically at 1 µL/min). The droplets were then flowing out of the collection vial towards the reinjection part of the chip where the fluorescence readout took place (typically at 10 -100 droplets per second). The delay time between the first droplet was produced and the first droplet was measured is typically 5 minutes.

Fluorescence measurement

The fluorescence signal of each droplet was measured as it was flowing through the microfluidic channel using a laser-induced fluorescence setup equipped with a 473 nm and a 532 nm laser (Drag-onLaser) and photomultiplier tube (PMT, Hammamatsu) in a system similar to the one described in Baret et al. 12 (Figure 2 (b,c)). The laser spot size was a few microns in diameter, much smaller than both the droplet and the channels (typically 100 microns). The optical filters and dichroics were purchased from AHF, the mounts and optical accessories from Thorlabs. Figure 3(a) shows a typical time trace of the droplets as they flow through the laser. For each droplet, we recorded the maximum fluorescence intensity in the green (fluorescein signal) and orange (resorufin signal) colors. We define the relative fluorescence unit as the measured voltage (U, from the PMT) rescaled by the Gain G as: RFU = U/G 6.8 . The exponent 6.8 depends on the PMT type, has been checked experimentally and corresponds to the data of the manufacturer. We recorded in realtime the droplet intensity as a function of time using a FPGA board (cRIO, National Instruments) controlled by a home-made Labview routine. The measurement of the maximum of the signal is optimized for fast processing with our FPGA board (up to several thousands of droplets per second processed in real time) but is sometimes sensitive to noise levels and laser instabilities leading to small variations in the measurements (of maximum 10 percent). The single point measurement provided a tool to measure the kinetics of equilibration of the emulsion.

Data processing

Data were post-processed with home-made scripts using standard toolbox in Matlab.

Experimental results

We performed the exchange experiments with resorufin. Two populations of monodisperse droplets This also indicates that the dye mainly remains in the aqueous compartments. We observed that the initial dye concentration did not significantly impact the time-scale of the process which justifies the rescaling of the concentration difference by the initial concentration c 0 (Figure 4(b)). In addition, the droplet size was not changing over the time-scale of the experiments and no transport of water was observed between the droplet: the transport of water is prevented by the high concentration of salts in the PBS buffer which was the same in all droplets 11,23 .

We performed the experiments for several surfactant concentrations above the CMC of the surfactant in HFE-7500. Figure 3(d) shows the time variation of the concentration difference over several hours of incubation. The rate of the equilibration increases with the surfactant concentration. All the experimental data collapse on a single master curve (Figure 4(a)) using a single time-scale τ e for each concentration indicating a similar process for the exchange at all concentrations above the CMC. τ e is inversely proportional to the surfactant concentration (Figure 4(b)) over almost two orders of magnitude in surfactant concentration. In contrast to other systems where the exchange is linked to bilayer formation 11 , our system shows the importance of the surfactant remaining in the continuous phase on the transport process. We observed no bilayers formed at any of these surfactant concentrations, therefore our results suggest that reverse micelles behave as nanocontainers increasing the solubility of the compound in the continuous phase 26 . The rescaled experimental data fit with an exponential decay only in the early time kinetics (Figure 4(a)). This exponential decay is not maintained at longer time-scale (closer to equilibrium), where a powerlaw behaviour is experimentally observed. In the region where the power-law is obtained, the fluorescence distributions of both population overlap. Notably, we can still measure the mean of each population due to the coding by fluorescein.

We also observed that the chemistry of the molecule influences the rate of exchange. As an example, fluorescein exchange is barely detactable over 24 hours. The exchange rates of sulforho-damine and dextran texas red are as well barely measurable on the time-scale of the experiments while rhodamine 101 is fully exchanged within minutes. This confirms previously published results 33 indicating that addition of sulfonate groups increases the chemical stability of the encapsulation. We also show that the time-scale τ e depends on the presence of Bovine Serum Albumin (BSA) in the droplet. The decrease of exchange by bovine serum albumin has been previously reported in mineral oils 28 . Here we measured here a 18 fold decrease of the exchange rate of resorufin by adding BSA in the droplet at a concentration of 25 mg/mL and the exchange rate is still proportional to the surfactant concentration (Figure 4 (b)).

One hypothesis, proposed previously in mineral oil emulsion, is that BSA acts as a barrier against exchange at the surface of the droplet 28 . We performed the following experiment to test this hypothesis and clarify the role of BSA in the control of exchange. All droplets contained the same resorufin concentration but the BSA (25 mg/mL) was present in only one type of droplet (Figure 4 (c,d

Modelling

The experimental results shown in figure 4(a) display a decay slower than exponential. For the quantity of interest, i.e. the difference between the mean values of the two subpopulations ∆c/c 0 , the system dimension affects the late kinetics with the power law ∆c/c 0 ∼ t -3/2 (Figure 4(a)). We developed a model to describe the exchange and explain the collective dynamics in the emulsion.

We consider now an ensemble of droplets separated by a porous membrane of thickness L 34 .

For simplicity, let us consider droplets arranged on the nodes of a three dimensionnal square lattice. The droplet on the node (k, l, m) contains molecules in a solvent at concentration c k,l,m and exchanges molecules with its 6 first neighbours only (4 neighbours in 2 dimensions, as displayed in Figure 5, 2 neighbours in 1 dimension). The area of contact between the compartments is s (a fraction of the droplet surface area S) and the volume of the compartment is V (Figure 5). The difference in resorufin chemical potential between the two compartments is responsible for the molecular flow with a flux J along each axis proportional to the concentration difference between neighbours δ c as described for membranes or bilayers 11,35,36 :

J = -Pδ c
J is expressed in mol/s/m 2 , the concentrations in mol/m 3 and the permeability P has the dimension of a speed (m/s). In this case the mass balance and Fick's law leads to a discrete version of the heat equation:

∂ c ∂ t = ∆c, (1)
where ∆ is a discrete Laplacian and t = t/τ is the dimensionless time with τ = V /(Ps). For periodic boundary conditions, the eigenvalues of the discrete Laplace operator are in three dimensions:

λ p,q,r = -4 × sin 2 π p N x + sin 2 π q N y + sin 2 π r N z
(for 0 ≤ p, q, r < N -1) and the eigenvectors

V p,q,r (k, l, m) = ω k p ω l q ω m r with ω p = exp(2iπ p/N x )
Due to the special form of the eigenvectors, the time evolution of the ensemble of N droplets is analytically expressed in the Fourier space as:

F (c) = F (c i )e λ p,q,r
t where c i is the initial distribution and F denotes the Discrete Fourier Transform. Since Eq. 1 is linear, we subtract in the following the mean value of concentration and normalize them to 1 to represent concentrations between -1 and 1, defining two subpopulations of droplets.

The smallest eigenvalue of the Laplace operator is -12 (in three dimensions) corresponding to the fastest decaying mode. The smallest eigenvalue is -8 in two dimensions, -4 for one dimension and -2 for a single droplet pair. It corresponds in all dimensions to an eigenmode of perfectly alternating array, where an empty droplet is surrounded solely by filled droplets and vice versa.

An alternated three dimensional packing would equilibrate following an exponential decay with a time scale τ/12. All other configurations decay more slowly.

The relaxation kinetics towards a uniform concentration depends on the excitation of the eigenmodes in the initial conditions. Experimentally, the mixing of droplets in the incubation vial leads to a random distribution of droplet positions. To investigate the effect of the disorder, we solved Eq. 1 for a completely random and uncorrelated initial condition. The scaling law can then be recovered analytically. Since the value of the initial concentrations is 1 or -1 depending on the subpopulation, ∆c/c 0 is proportional to the scalar product of c by c i . Using Parseval's theorem, this scalar product can be computed in the Fourier space. In 3-d, we have:

∆c c 0 ∝ ∑ p,q,r
|F (n 0 )| 2 e λ p,q,r t Due to the completely random initial condition, all the modes are statistically excited with the same weight (which is an initial condition) and ∆c/c 0 ∝ ∑ p,q,r e λ p,q,r t . At long time, only the modes with the smaller eigenvalues remain, so that

λ p,q,r ∼ -4π 2 ((p/N x ) 2 + (q/N y ) 2 + (r/N z ) 2 )
and the sum can be extended to +∞. For large systems, we can use a continuum approximation, and

∆c c 0 ∼ +∞ 0 +∞ 0 +∞ 0 e -4[(π p/N x ) 2 +(πq/N y ) 2 +(πr/N z ) 2 ]t dpdqdr
which scales as t -3/2 (in 3d). Figure 5(e) shows a summary of the model data in the three dimensional case and the corresponding asymptotic regimes. Our model can also be solved in 1-d and 2-d cases and for any set of initial conditions using similar schemes. Similarly, we also recover analytically that ∆c/c 0 scales in 1-d like as t-1/2 , and in 2-d as t-1 .

At short-times, an argument similar to the one at large time-scales leads to ∆c/c 0 ∼ ∑ p,q,r e λ p,q,r t .

Expanding the exponential close to 0 leads to ∆c/c 0 ∼ 1λ t, where λ is the average value of the eigenvalues, which is equal to -6 in three dimensions (-4 in two dimensions and -2 in one dimension). Therefore the early time kinetics in a random initial distribution is an exponential decay two times slower than the exponential decay of the fastest mode (the alternating case) of the corresponding dimension.

In summary, the relaxation of the chemical potential at the emulsion level is a function of the initial distribution of the droplet in the packing. The maximum rate is obtained for an alternating packing and corresponds to the rate of exchange at the droplet pair level (corrected by the number of nearest neighbours). The early time kinetics of a random distribution is twice slower but has the same order of magnitude, while alternatively, the rate of exchange of an emulsion made of two layers would equilibrate with a much longer time-scale, that can be several orders of magnitude longer. These results are generic and do not rely on the exact description of the exchange mechanism at the droplet level.

Discussion

Our experiments show how the packing and the random organisation of the droplets affect the kinetics of equilibration of the emulsion. A single time-scale controls the exchange and we measured how this time-scale is influenced by surfactant concentration, the chemistry of the molecules to be transported and presence of additives such as BSA in the dispersed phase. We want here to relate our experiments to models of membrane permeability 23,35,36 .

When the rate-limiting step is diffusion in the membrane (case I of Zwolinksi et al. 35 ), the membrane permeability is expressed as P = KD/L 35 , K being the partition coefficient of the dye between the reservoir and the membrane, D the diffusion coefficient of the dye and L the membrane thickness. One of the control parameter is therefore the partition coefficient between the compartment and the membrane. Non-fluorinated organic molecules, such as resorufin, are typically highly insoluble in fluorinated oils similar to those used for the continuous phase in this study 29 . We have not been able to directly measure the partition coefficient between the aqueous phase and the oil, or the oil-surfactant mixtures as, (i) the optical properties of the dye seem to differ in water and in oil and (ii) the presence of surfactant makes such bulk measurement technically challenging. Furthermore, performing an experiment with BSA in half of the droplets provided a means to indirectly address the role of the partition coefficient. Starting from droplets which all have an identical concentration of resorufin, but where half of the droplets contained 25 mg/mL of BSA, the concentration of resorufin in the two populations diverged until, at equilibrium, there was 13 no BSA. Therefore adding BSA provides a means to modify the parameter K without modifying the dye or the continuous phase. Exchange experiments performed with 25 mg/mL BSA in all droplets and with resorufin in only half led to a decrease of exchange rate by a factor of about 18 which matches quantitatively the corresponding decrease of K. Our experiments are therefore consistent with a diffusion-limited permeability model. This mechanism is also consistent with the results obtained with other dyes showing a decrease of exchange when sulfonated groups are added 33 : if the solubility in the aqueous phase increases, the rate of exchange decreases.

The decrease of solute exchange mediated by BSA has here been quantitatively measured with fluorinated oils, but was previously reported for non fluorinated oils 28 . One of the initial explanation for the effect of BSA is the creation of a protein layer at the inner surface of the droplets which provides a barrier in a capsule-like manner as suggested earlier for mineral oils 28 . This layer would decrease the permeability of the medium between the droplets 37 . In this case, the ratelimiting step would be the diffusion through the BSA layer (case II or III of Zwolinksi et al. 35 ).

However with such picture one would not expect that (i) the surfactant concentration plays a role in the presence of BSA and that (ii) starting from equilibrated concentrations, the final state shows a concentration difference. Our experiments show unambiguously that BSA acts on the solubility of the dye (a thermodynamic property) rather than on kinetic parameters such as energy barriers.

The quantitative link between the equilibrium concentration ratio of resorufin in aqueous droplet with and witout BSA and the exchange rate with BSA is consistent with a model of porous membranes where the permeability is controlled by partition coefficient between the aqueous phase (the droplets) and the continuous phase (fluorinated oil containing surfactant). We speculate that the change in partition coefficient of resorufin caused by BSA may be due to the remarkable capacity of serum albumins to bind a structurally diverse set of small molecules, including many drugs [START_REF]All about albumin: biochemistry, genetics and medical applications[END_REF] . The binding of resorufin, or other solutes, to BSA could potentially increase solubility in the aquous phase and/or decrease solubility in the oil/surfactant phase.

Recently the permeability of a bilayer membrane has been linked to its fluidity 11 . In our case, the surfactant does not form bilayers and the micelle concentration in the oil phase plays a significant role in the transport, the exchange rate being proportionnal to micelle concentration. It indirectly shows that the solubility in the continuous phase is increasing linearly with surfactant concentration. Therefore the likely mechanism for transport is through the solubilisation of the dye molecule into the reverse micelles (above the CMC, the free surfactant concentration is constant) which can act as relatively large (∼ 100 nm diameter) carrier of molecules between droplets. The shape and size of the micelles could also be an important control parameter in the transport and further experiments investigating the structure of these micelles would be of interest.

Finally, for a given emulsification system, the time-scale τ = V /(Ps) controls the exchange, with a prefactor corresponding to the number of nearest neighbours. We anticipate that the area controlling the transport is proportional to the droplet area R 2 . However, in the case of bilayers, there is no reason to assume that P also scales with the droplet size. Therefore τ probably has a linear dependence with the droplet size τ ∼ R/P and we expect that a decrease of droplet size will increase the rate of exchange. On the contrary, when no bilayers are formed P integrates a geometric parameter such as the droplet-droplet distance and the scaling will be different. Further experiments investigating these scalings will be of interest. The link between the exchange rates in the different dimensions is determined by the physical details of the membrane between droplets.

Microfluidic experiments provide a means to control droplets to study isolated pairs or generate controlled packing in 1d or 2d and 3d. By measuring the exchange kinetics in the different dimensions and varying the number of nearest neighbours, it will be possible to extract microscopic information on the membrane (for example what is the exact description of P and s).

Microfluidic systems have a huge potential to address the remaining open questions by the possible on-chip control of packing geometries, droplet spacing or sizes for a deeper understanding of such transport phenomenae. We plan as a next step to study the effect of the packing geometry on the dynamics of the relaxation as well as the effect of targeted or random droplet rearrangements: Inducing re-ordering of droplets while diffusion occurs might lead to interesting kinetics -especially when the droplet reshuffling can be controlled [START_REF] Surenjav | [END_REF] . Microfluidics will provide, in the future, a very powerful platform for the study of such transport problems by the possible control of droplet packing and lateral ordering of the droplets [START_REF] Surenjav | [END_REF][40][41] . Finally, the mechanisms involved in this study are generic (a diffusive process) and could potentially be applied to a better control of collective effects in droplet population, for example, in active emulsions [42][43][44] , to control signaling in cell populations (chemotaxis, quorum sensing) entrapped in different ecological niche (here the droplets), or to design new systems for screening compounds by affinity. The ability to control the transport by additives selectively encapsulated in specific droplets in the emulsion would also open the door to applications for separation.

Conclusions

In summary, microfluidics provides the tools to produce calibrated emulsions and to analyse millions of droplets at very high-throughput. We have measured experimentally the equilibration of the chemical potential of chemicals in an emulsion initially composed of two sub-populations of different dye concentrations. We measured the exchange timescale for a three dimensional packing of randomly mixed droplets. We observed that the exchange timescale in a random packing is controlled by the number of reverse micelles in the continuous phase and relates to a classical description of porous membranes through partition coefficient between the dispersed and continuous phase. We showed that -due to the disorder in the initial condition -the macroscopic kinetics are not a simple exponential decay but a power-law and that the initial conditions have a major influence on the behavior of the system. We also showed how the initial distribution of droplets can affect the time-scale of the exchange. The early time kinetics in a mixed emulsion is only twice slower than the maximum rate, obtained for a perfectly organized alternating packing. Most of the important physics of the relaxation is contained in our simple model: the equilibration process is a classical diffusion process through a Laplace equation from which emerge a dimension-dependent power-law. The transport -as a passive process -leads to equilibration of concentration and can be modulated by additives such as BSA. Selective transport is obtained by controlling the BSA en-capsulation in some droplets to concentrate the compound of interest, leading to a versatile control of chemicals in emulsions. The background count rate of the pure oil is given by the horizontal full line. The grey area represents the standard deviation for 10 measurements. At low surfactant concentrations (below 0.02%, the scattered count rate is in the oil background). Above 0.02%, the count rate of the diffused light is quasi-linear (exponent 1.13). Inset: The radius ρ of reverse micelles measured by correlation of the scattering signal is slightly increasing around ρ = 80 nm up to 120 nm at 5%. The first concentration for which a correlation is obtained is 0.02% although in some runs at 0.02% the background is too large to obtain a proper correlation. The diamond has been obtained for a 1 µM resorufin initial concentration. The open circle correspond to the rescaling in presence of Bovine Serum Albumine (the dashed line and the dotted-dashed line correspond to a linear relationship and the ratio of both is 18). (c) Exchange kinetic at 2% surfactant concentration when all droplet initially contain 5 µM resorufin but only half contain 25 mg/mL BSA. The equilibrium shows a 15 fold ratio in concentration of resorufin between both types of droplets, c w is the resorufin concentration in the droplet with BSA and c wo without BSA (d). This ratio closely matches the ratio of 18 observed in the exchange rate of (b).

  From the raw data, we generated the time distribution of the fluorescence by binning typically 2000 droplets (Figure 3(b)). We generated color plots corresponding to the time evolution of the histogram (Figure 3(c)). From the signals, we first analysed the green fluorescence to distinguish the two populations of droplets. We checked that fluorescein signals were constant at the time-scale of our experiments. The orange fluorescence of both population was then fitted by Gaussian distributions to obtain the means m 1 and m 2 of each droplet population. In the linear regime of fluorescence (as it is the case here), the fluorescence was proportional to dye concentration. Therefore the value m 2m 1 is proportional to the difference of concentration of the dye. By rescaling the data by the initial value, we obtained a measurement of the concentration difference rescaled by the initial concentration: ∆c/c 0 = (m 2m 1 )/(m 2m 1 ) 0 . This parameter -equal to 1 at the beginning of the experiment and to 0 when the concentrations are equilibrated -provided means to compare experiments performed for different dyes and different initial ratio of dye concentrations.

  (0.3 nL) with different concentration of resorufin encoded by fluorescein are mixed. We measured the kinetics of equilibration of the chemical potential of the resorufin by fluorescence. Figure 3(b) and (c) show the typical time distribution of the droplet fluorescence over a few hours for C = 0.5 % w/w . Starting from two distinct populations of droplets, the equilibrium was reached after several hours. The average value of the fluorescence signal of the two populations was constant throughout the measurement indicating that the dye is preserved (no chemical modification nor bleaching).

  )), those droplets being encoded with a 150 nM fluorescein code (droplet containing no BSA also contain 50 nM fluorescein). We measured the variation of the resorufin fluorescence signal as a function of time. In this experiment, we also observed the exchange of resorufin leading to an a priori counter-intuitive result: the difference in fluorescence is increasing during incubation (with a time-scale of order 20 minutes). At equilibrium, the fluorescence signals of the two types of droplets are different by a factor of 15 (Figure4 (d)) and the most fluorescent droplets are those containing the BSA. The concentration of resorufin in the droplets containing BSA is 15 fold larger than in droplets containing no BSA. Our experimental results clearly contradicts the assumption that BSA creates a thin viscoelastic layer at the interface of the droplets acting as a kinetic barrier. If the BSA only acts as a barrier to exchange, the chemical potential of the resorufin in the droplets being initially at equilibrium, the fluorescence signal should remain constant. The time-scale of the exchange process is similar to the time-scale of exchange previously measured for the corresponding surfactant concentration (2%) which indicates that the mechanism of transport by surfactant is not affected by the presence of BSA. Finally, the ratio between the equilibrium fluorescence of the two droplet populations matches quantitatively the decrease of exchange rate obtained in Figure4(b). In this particular case, BSA is efficiently used to extract and concentrate resorufin.

Figure 1

 1 Figure 1 Surfactant solutions (a) The surfactant used in this study is a Krytox-Jeffamine block-copolymer (x + z ≈ 3.6, y ≈ 9.0, n =35-40) solubilized in fluorinated oil (HFE-7500) (b). The fluorinated oil is used for its low ability to solubilize organic compounds. (c) Equilibrium surface tension for the water/HFE interface in presence of the surfactant at different concentrations. We checked the reproducibility of the synthesis with two batches of surfactant (+ and ×). (d) Dynamic Light Scattering on the oil -surfactant solution.The background count rate of the pure oil is given by the horizontal full line. The grey area represents the standard deviation for 10 measurements. At low surfactant concentrations (below 0.02%, the scattered count rate is in the oil background). Above 0.02%, the count rate of the diffused light is quasi-linear (exponent 1.13). Inset: The radius ρ of reverse micelles measured by correlation of the scattering signal is slightly increasing around ρ = 80 nm up to 120 nm at 5%. The first concentration for which a correlation is obtained is 0.02% although in some runs at 0.02% the background is too large to obtain a proper correlation. Combining both methods leads to a CMC measurement in the range to 0.02 -0.04% w/w .

Figure 2

 2 Figure 1 Surfactant solutions (a) The surfactant used in this study is a Krytox-Jeffamine block-copolymer (x + z ≈ 3.6, y ≈ 9.0, n =35-40) solubilized in fluorinated oil (HFE-7500) (b). The fluorinated oil is used for its low ability to solubilize organic compounds. (c) Equilibrium surface tension for the water/HFE interface in presence of the surfactant at different concentrations. We checked the reproducibility of the synthesis with two batches of surfactant (+ and ×). (d) Dynamic Light Scattering on the oil -surfactant solution.The background count rate of the pure oil is given by the horizontal full line. The grey area represents the standard deviation for 10 measurements. At low surfactant concentrations (below 0.02%, the scattered count rate is in the oil background). Above 0.02%, the count rate of the diffused light is quasi-linear (exponent 1.13). Inset: The radius ρ of reverse micelles measured by correlation of the scattering signal is slightly increasing around ρ = 80 nm up to 120 nm at 5%. The first concentration for which a correlation is obtained is 0.02% although in some runs at 0.02% the background is too large to obtain a proper correlation. Combining both methods leads to a CMC measurement in the range to 0.02 -0.04% w/w .

Figure 3 Figure 4

 34 Figure 3 Fluorescence measurement and Analysis (a) Signals measured on the PMT at droplet production. Droplets containing 10µM of resorufin (droplets 1,2,5,7 and 8) are encoded with a high concentration of fluorescein (150 nM). The droplet initially containing no resorufin (3,4 and 6) are encoded with a low fluorescein concentration (50 nM). (b) Histogram of resorufin fluorescence at different reinjection time for C = 0.5% w/w of KryJeffa(1) . The histograms are obtained from about 2000 droplets each. (c) Full time histogram (colorplot) of the orange fluorescence recorded over ∼ 8 hours for a mixed emulsion of resorufin in HFE-7500 with C = 0.5% w/w of KryJeffa(1) -surfactant showing the equilibration of resorufin chemical potential. (d) Measurement of the concentration difference of the two populations of droplets for increasing surfactant concentration. The kinetics are faster for larger surfactant concentrations.

Figure 5

 5 Figure 5 Model for molecular exchange. (a) For a porous membrane separating two compartments, the flux J depends on the concentration difference and geometry parameters (s, L). For droplet pairs we use a similar description (b), the flux being proportional to the concentration difference. In an ensemble of droplets, the exchange occurs with the nearest neighbors: (c) corresponds to a one-dimensional case, (d) the two-dimensional case on a square lattice. (e) The model is solved for a three-dimensional packing on a square lattice. Results of the model and asymptotic regimes in the three-dimensional case.

Aknowledgments

The authors acknowledge the financial support of Region Alsace, of the 'Fondation d'Entreprise EADS' and of the European Commission Framework Program 6 (MiFem). Financial support from the Max-Planck Gesellschaft and from the Deutsche Forschunggemeinschaft through the Sonderforschungbereich 755 (Nanoscale Photonic Imaging) are also acknowledged. The authors are grateful to C. Rick, C. Ben Salem and Z. El Harrak for help and assistance as well as B. Maury, S. Herminghaus, J. Bibette and N. Bremond for helpful discussions.