Vincenzo Di Paola
email: vincenzo.dipaola@edu.unige.it

Alexandre Goldszstejn

Matteo Zoppi
email: matteo@unige.it

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

Alexandre Goldsztejn
email: alexandre.goldsztejn@ls2n.fr

Analytic Center Based Tension Distribution for Cable-Driven Platforms (CDPs)

Keywords: Parallel Robots, Cables, Actuation Redundancy, Tension Distribution, Robustness, Continuous

A redundant Cable-Driven Platform (CDP) is composed of m cables that exceed the Degree of Freedom (DOF) of the end-effector. The choice of tension along the cables admits infinite solutions. This paper proposes the use of the Analytic Centre to solve the tension distribution problem. Adopting this technique allows finding tensions far from the tension limits namely, robust as well as tension profiles continuous and differentiable in time. The continuity, differentiability and uniqueness of the solution is also proved. Moreover, the possibility of including non-linear constraints acting on the tensions (e.g. friction) is a further contribution. The computational time of proposed approach is compared to the existing techniques to asses its real-time applicability. Finally, several simulations using several CDPRs' architectures are reported to demonstrate the method's capabilities.

I. INTRODUCTION

C ABLE-Driven Platform (CDP) constitutes an interesting class of cable-driven robots, which have experienced a massive growth in the last decade. Indeed, within this class of robots one can find the widely established Cable-Driven Parallel Robots [START_REF] Bruckmann | Cable-driven parallel robots[END_REF] (CDPRs) and the newcomers Aerial Cable Towed Systems [START_REF] Villa | A survey on load transportation using multirotor uavs[END_REF] (ACTSs) or a hybrid version between CDPRs and ACTSs also known as Hybrid-ACTSs [START_REF] Jamshidifar | Static workspace optimization of aerial cable towed robots with land-fixed winches[END_REF]. The rapid development of CDPs is due to their inherent capabilities such as high tracking speed, large workspace, reconfigurability and modularity, which make them suitable for a vast category of applications such as crane-like applications [START_REF] Albus | The nist robocrane[END_REF], [START_REF] Merlet | A portable, modular parallel wire crane for rescue operations[END_REF], rehabilitation [START_REF] Rosati | Design, implementation and clinical tests of a wire-based robot for neurorehabilitation[END_REF], [START_REF] Surdilovic | String-man: Wire-robot technology for safe, flexible and human-friendly gait rehabilitation[END_REF], filming [START_REF] Cone | Skycam: An aerial robotic camera system[END_REF], aerial load transportation to cover long distances [START_REF] Michael | Cooperative manipulation and transportation with aerial robots[END_REF], post-disaster scenarios [START_REF] Bosscher | A concept for rapidly-deployable cable robot search and rescue systems[END_REF] and so on.

Generally speaking, the presence of cables brings advantages such as reduced inertia and low cost while making them collaborative with humans. Despite mentioned advantages, there is an added complexity in operating with cables. Indeed, guaranteeing cable tautness and smooth variations of the tension values during the execution of a task poses several problems. The management of cable tensions is still a topic of ongoing research. When the Degree of Redundacy DoR1 ≥ 1 namely redundant CDPs are considered, infinitely many solutions to the tension distribution problem may exist. In this view, one chance of solving this problem is to formulate an optimisation problem and solve it iteratively.

The 1-norm of the tension vector τ is the protagonist of several previous works [START_REF] Borgstrom | Rapid computation of optimally safe tension distributions for parallel cable-driven robots[END_REF]- [START_REF] Oh | Cable-suspended planar parallel robots with redundant cables: controllers with positive cable tensions[END_REF]. Linear Program (LP) are notable for their fast convergence, which makes them suitable for real-time applications. Though, the optimal point is always placed at the edge of the tension box and, between two iterations, it can jump from one edge of the polyhedron to the other resulting in discontinuities between two successive solutions. The same can happen also with the ∞-norm [START_REF] Gosselin | On the determination of the force distribution in overconstrained cable-driven parallel mechanisms[END_REF] thus, to overcome this issue and guarantee continuity between solutions, the p-norm 1 < p < ∞ can be used [START_REF] Verhoeven | Tension distribution in tendon-based stewart platforms[END_REF]. In fact, practically, a natural extension consists in using the 2-norm of the tension vector [START_REF] Taghirad | An analytic-iterative redundancy resolution scheme for cable-driven redundant parallel manipulators[END_REF], [START_REF] Agahi | Redundancy resolution of wire-actuated parallel manipulators[END_REF] and in solving a Quadratic Program (QP) whose accessibility makes it popular. Even so, other solutions were explored for example, in [START_REF] Hassan | Optimization of actuator forces in cablebased parallel manipulators using convex analysis[END_REF], [START_REF]Analysis of bounded cable tensions in cable-actuated parallel manipulators[END_REF] the Dykstra alternating-projection algorithm is used to solve the minimum 2-norm tension distribution. Another alternative to tackle the tension distribution problem, ensuring continuity, takes the barycenter of the polyhedron as the optimal solution [START_REF] Mikelsons | A real-time capable force calculation algorithm for redundant tendon-based parallel manipulators[END_REF]. Evidently, the criteria are different and their choice depends on the assigned task. Indeed, for example, the 2norm methods can guarantee reducing power consumption whereas barycenter is said to be safe since the solution is far from cable tension limits. Besides, the Improved-Closed Form method [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF], built on its predecessor [START_REF] Pott | Closed-form force distribution for parallel wire robots[END_REF], tries to reduce the method complexity by maintaining real-time efficiency, generality w.r.t. the DoR and continuity of the tension profiles. Even so, these methods do not guarantee the convergence in the entire Wrench Feasible Workspace (WFW) [START_REF] Bouchard | On the ability of a cabledriven robot to generate a prescribed set of wrenches[END_REF]. In addition, since the Improved-Closed Form approach set outbound cable tensions to their maximum or minimum value, the smoothness of tension profiles can be compromised. Analogously, with the intention of maximising the applicability of the method within the WFW without degrading continuity and algorithm speed, in [START_REF] Katharina Müller | Analysis of a real-time capable cable force computation method[END_REF] they take their inspiration from [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF] and developed the Improved Puncture Method. A strategy for tension distribution dealing with CDPRs operating beyond their WFW is also explored in [START_REF] Côté | A tension distribution algorithm for cable-driven parallel robots operating beyond their wrench-feasible workspace[END_REF]. Actually, many works were proposed to establish an opportune set of tensions to be used during the fulfilment of a given task. Anyway, only one [START_REF] Gouttefarde | A versatile tension distribution algorithm for n-dof parallel robots driven by n + 2 cables[END_REF] introduces a versatile tension distribution algorithm, which allows computing several types of mentioned techniques. Nevertheless, the introduced algorithm can be applied only for CDPRs with, at most, DoR = 2.

In this context, this paper aims to introduce the Analytic Centre method. This new technique provides robust (i.e. safe) continuous and differentiable tension profiles for any DoR. Another contribution consists in the ability to incorporate time-varying non-linear constraints acting on the tensions. This latter increases even more its range of applicability and generality to other types of CDP as MCDPRs [START_REF] Tahir Rasheed | Tension distribution algorithm for planar mobile cable-driven parallel robots[END_REF], [START_REF] Rasheed | Wrench-feasible workspace of mobile cable-driven parallel robots[END_REF] and moreover, it constitutes another step towards the subject of Physical-human Interaction (PhI). Indeed, the possibility of limiting the tensions of specific cables over time (i.e. during specific parts of the trajectory) materialises the possibility of having the human close to certain cables and during specific phases of a given task.

This paper is structured as follows: in Section II the CDP model and the main definitions necessary to state the tension distribution problem are recalled. A motivating example opens Section III followed by the definition of the Analytic Centre in Section IV. The continuity and differentiability are rigorously proved as well as the existence of a unique solution for the parametric optimization problem. Subsequently, Section V includes simulations with various CDPR architectures while last Section VI draws the conclusion and outlines future works.

II. CDP MODEL

In this section, the main equations and symbols necessary to describe the CDP are hereby reported.

The static or dynamic equilibrium of a platform in the space, guided by m cables is governed by the following equation

Wτ + w e = 0, (1)
where, in general, w e ∈ R 6 is the external wrench which also takes into dynamical actions applied to the platform, τ ∈ R m is the cable tensions vector and the terms W ∈ R 6×m is the wrench matrix which is defined as

W = u 1 . . . u m b 1 × u 1 . . . b m × u m , (2)
here u i ∈ R 3 represents the ith cable direction unitary vector and b i ∈ R 3 represents the ith attachment point on the platform. Thus, if DoR ≥ 1, there exist infinite solutions of Eq.( 1) grouped in the following set

Σ = τ | Wτ + w e = 0 . (3)
However, to maintain the equilibrium of the platform the cable tension limits have to be taken into account. Hence, the m-dimensional convex hypercube Π that defines the domain of the feasible tensions is

Π = τ | 0 < τ ≤ τ ≤ τ , (4)
where τ , τ ∈ R m,+ are positive tension vectors limits containing the ith lower and upper cable tension limits as a components which, without loosing generality, will be considered equal to each other, respectively. Consequently, the set of feasible solutions Γ satisfying both Eq.( 1) and Eq.( 4), is

Γ = Σ ∩ Π. (5)
Practically, the cable tensions can be computed as

τ = -W † w e + Nλ = τ p + τ g , (6)
where W † = W T (WW T ) -1 is the Moore-Penrose pseudoinverse matrix of W, N ∈ R m×(m-n) contains the vectors that span the kernel of W, λ ∈ R (m-n) belongs to the polytope Λ Eq.( 7) and τ p , τ g are the particular and general solution of Eq.( 1), respectively. For sake of clarity the definition of the convex polytope Λ is reported in

Λ := {λ ∈ R (m-n) | τ min ≤ -W † w e + Nλ ≤ τ max }. ( 7
)
It turns out that to cope with the cable intrinsic property, guaranteeing the controllability of the platform during the tasks, a natural way to deal with Eq.( 6) consists in solving an optimization problem. So far, several approaches were presented in the literature. Therefore, a motivation and an overview of the context in which this contribution lies is given in the next section.

III. MOTIVATING EXAMPLE

A. Description

This section consider an academic example that consists in distributing two tensions τ = (τ 1 , τ 2 ) belonging to a box Π = τ | 0 < 10 ≤ τ ≤ 100 with corners τ = (10, 10) and τ = (100, 100). The level-sets of different tension distribution optimization problems are displayed in Figure 1. Tensions are required to satisfy the static equilibrium constraint Wτ + w e = 0 with W ∈ R 1×2 and w e ∈ R 1 . Three test cases are considered: the first takes W = (-7, 20) and w e = (-1790) N and it is displayed in black in Figure 1. The second takes W = (-1, 50) and w e = (-945) N and it is displayed in blue in Figure 1. The third takes W = (1, 50) and w e = (-1055) N and it is displayed in green in Figure 1. The last two sets of data are very similar to each other, being representative of close poses or of the effect of uncertainties on a pose, and are presented to illustrate the sensitivity of the solution given by Tension Distribution Algorithms (TDAs).

In all the cases, the barycenter of the feasible polytope Λ coincides with the midpoint of the two intersections between the line Wτ + w e = 0 and the borders of Π.

B. State-of-the-art optimization based TDAs

TDAs based on minimizing the 2-norm have been widely investigated and used. They rely on solving the linearly constrained quadratic problems:

argmin Wτ +we=0 τ ∈Π ∥τ ∥ 2 and argmin Wτ +we=0 ∥τ -τ ∥ 2 (8) (a) (b) (c) (f) (e) (d)
Fig. 1. Level-sets: (a) the 2-norm with respect to the origin (0, 0), (b) the 2-norm with respect to the center of Π, (c) the 5-norm within Π, (d) the ∞-norm within Π, (e) the distance to the closest boundary of Π and (f) the logarithmic barrier over Π. Three test cases are displayed in black, blue and green respectively. From Lagrange optimality conditions, if a minimizer lies strictly inside the cable tension box then the level-set of the cost function and the linear subspace Wτ + we = 0 are tangent at this minimizer. Using generalized gradients, this tangency condition is satisfied at corners of non-smooth level-sets in the lower left and middle graphics. with τ = 1 2 (τ + τ ) the midpoint of Π. The first allows minimizing the required energy to the actuators, hence it offers a limited stiffness. Furthermore, its minimizer often lies on the boundary of Π and is therefore not considered in the context of robust tension distribution. The level-sets of the cost function ∥τ ∥ 2 and the corresponding minimizers are shown in Figure 1-(a). It is possible to see that green and blue tension vectors are close to each other, indicating a low sensitivity of the investigated TDA. Minimizing ∥τ -τ ∥ 2 is meant to compute tension vectors as close as possible to the center of Π. It is therefore attracting for pursuing robust tension distribution. Using the 2-norm allows for a Closed Form solution of Eq.( 8), see [START_REF] Pott | Closed-form force distribution for parallel wire robots[END_REF]. This Closed Form solution is additionally attractive because its computation requires only to evaluate the pseudoinverse of W. The green and blue solutions displayed in Figure 1-(b) are seen to be very close to each other, showing again a low sensitivity of this TDA.

Since the inequality constraints are not taken into account in the optimization problem, its solution can be outside Π while some feasible solution exists, and indeed the black solution is not visible in Figure 1-(b) because it is outside Π (see also Figure 2 in [START_REF] Pott | Closed-form force distribution for parallel wire robots[END_REF] for a more detailed discussion). This results in a reduced WFW coverage. A tentative post-process aiming setting to the lower or upper limit some tensions computed by the Closed Form formula that exceed the limits was introduced in [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF], but with no sound of theoretical foundation for the convergence and continuity. Furthermore, this corrected Closed Form tension vector is meant to lie on the boundary of Π, hence not satisfactory with respect to robustness.

Remark III.1. Correcting the Closed Form formula (or directly) solving Eq.( 8) with the inequality constraints τ ∈ Π using a numerical solver has not been considered in the literature yet to the best of authors knowledge. Its numerical optimization is as complex as with τ = 0, and is therefore also attractive. However, when the Closed Form solution is outside Π, the corrected constrained solution will lie on Π boundary, hence it is not satisfactory in term of robustness.

Remark III.2. Using p-norms, including p = ∞, was discussed in [START_REF] Tobias Bruckmann | Calculating force distributions for redundantlyactuated tendon-based stewart platforms[END_REF] while the case p = 4 was extensively discussed in [START_REF] Gosselin | On the determination of the force distribution in overconstrained cable-driven parallel mechanisms[END_REF] where a Closed Form solution is also given. The level-sets of the ∞-norm with respect to the center of Π are shown on Figure 1-(d). The tangency condition (i.e. optimality) is satisfied when the linear subspace meets a corner of the level-set, as illustrated by the three solutions represented. Because the level-sets are straight lines, cases where the linear subspace Wτ + w e = 0 is parallel to the level-sets have an entire segment of solutions, and an infinitesimal change in the data of the problem will cause an abrupt transition to one or the other segment endpoints. This occurrences is illustrated on Figure 1-(d) by the high sensitivity of the blue and green neighbor problems whose solutions are characterized by a big distance. Intermediate behaviors are obtained considering pnorms with 2 < p < ∞: Figure 1-(c) shows the level-sets of the 5-norm. The black solution computed is now inside Π.

There is no straight level sets anymore, but the presence of almost straight level-sets show a strong sensitivity as well, as illustrated by the large distance between the blue and green solutions.

Aiming maximizing the robustness of the TDA, maximizing the distance to the closest face of Π was proposed in [START_REF] Borgstrom | Rapid computation of optimally safe tension distributions for parallel cable-driven robots[END_REF]. This leads to the following optimization problem:

argmax Wτ +we=0 τ ∈Π d(τ ), (9)
where d(τ ) = min{τ -τ 1 , τ 1 -τ , τ -τ 2 , τ 2 -τ } is the distance to the closest face of Π. This optimization problem is classically reformulated to an easy-to-solve Linear Program (LP). Figure 1-(e) shows the level-set of d(τ ). In spite of its advantage in term of robustness, this approach shows a very high sensitivity illustrated by the large distance between the blue and green solutions. This happens when the linear subspace Wτ + w e = 0 is close to be parallel to an axis. Although in practice the exact time where the linear subspace is exactly parallel to an axis is in general not at a timestep where the optimization problem is solved, this leads to potentially discontinuous tension profiles as detailed in [START_REF] Borgstrom | Rapid computation of optimally safe tension distributions for parallel cable-driven robots[END_REF].

Remark III.3. Figures 1-(d)(e)
show a surprising coincidence between the level-sets of the ∞-norm with respect to the center and the distance to the closest face. This coincidence leads to the same minimal and maximal tension vectors. As far as the authors know, this coincidence was not made explicit in this context, and actually holds in generals when tension limits are independent of cables2 , i.e., τ i = τ j = τ and τ i = τ j = τ . The distance to the closest face of Π is then

d(τ ) = min{τ -τ 1 , τ 1 -τ , . . . , τ -τ m , τ m -τ }. ( 10
)
In the following, the constraints Wτ + w e = 0 and τ ≤ τ ≤ τ are implicit in the optimization problems to lighten the notation. Recall that argmax f = argmin -f and that adding the radius r = 1 2 (τ -τ ) to the cost function does not change its argument. Furthermore, the radius r and center τ = 1 2 (τ + τ ) are related by τ -r = τ and τ + r = τ . Then basic manipulations show that:

argmax d(τ ) = argmin -d(τ ) + r (11a) = argmin max{τ 1 -τ , τ -τ 1 , . . . , τ m -τ , τ -τ m } (11b) = argmin max{|τ 1 -τ |, . . . , |τ m -τ |} (11c) = argmin ∥τ -τ ∥ ∞ . (11d)

C. The Analytic Centre

A great advantage of the TDA corresponding to the distance to center with 2-norm and no inequality constraint is that the absence of inequality constraints allows an efficient resolution. In this special case a Closed Form solution is available but, in general, numerical optimization is easier and more efficient when only linear equality constraints are enforced: basically, inequality constraints require determining which inequality is active at minimizer, entailing more complex algorithms. One escamotage in optimization foresee the usage of logarithmic barrier functions, which change inequality constrained minimization problems of f (x), subject to inequality constraints g i (x) ≤ 0, to unconstrained minimization problems consisting in minimizing f (x) -c i log -g i (x) . In the context of inequality constrained optimization, the parameter c is meant to converge to zero during the resolution in order to reduce the impact of the barrier function on the objective, see Figure 2.

Remark III.4. The logarithmic barrier is preferred to other barrier functions like 1 g(x) or 1 g(x) 2 because it enjoys the technical property of being self-concordant, a theory introduced in the late 80's by Nesterov and Nemirovski [START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF] which allows bounding the number of Newton iterations required to reach a prescribed accuracy on the minimum. Such a good convergence characterization of the numerical optimization is of critical importance for real time usage of the TDA.

For a fixed value of the barrier coefficient c, the barrier actually attracts the solution inside the solution set. Hence, one aim of this paper is to use this property to enforce the robustness on the tension solution. This can be achieved considering the so-called Analytic Centre associated with the constraints τ ∈ Π and Wτ + w e = 0, which is defined as:

argmin Wτ +we=0 ϕ(τ ) (12)
with

ϕ(τ ) = 2 i=1 -log(τ -τ i ) -log(τ i -τ ). (13)
The Analytic Centre definition, as well as many properties and related methods of numerical optimization, can be found in [START_REF] Boyd | Convex optimization[END_REF]. As seen in Eq.( 12), logarithm barrier function enforce τ to be strictly inside Π.

Remark III.5. In the framework of convex optimization, it is classical to consider the constraint τ ∈ int(Π) implicit in Eq.( 12) with ϕ(τ ) = +∞ when logarithms are not defined. In particular, this emphasizes that optimization methods dedicated to problems with linear equality constraints and no inequality constraints are used to solve Eq.( 12).

The level-sets of ϕ(τ ) are shown in Figure 1-(f). One can observe that the Analytic Centre optimization problem offers a much better compromise between the 2-norm to the center and the ∞-norm to the center of Π3 than the 5-norm does: first, contrarily to the 2-norm to the center and similarly to the ∞-distance to the center of Π, the Analytic Centre is inside Π under the weak assumption that the linear subspace Wτ + w e = 0 intersects the interior of Π4 , see the black solution for a typical case where the minimal 2-norm is outside Π while the analytic center is inside. Second, the level sets are similar to circles in a neighborhood of the center of Π, hence showing a low sensitivity in this area, as illustrated by the blue and green solutions in Figure 1-(f), which can be compared .5 -0.5 0 with Figure 1-(c) where the blue and green solutions are far away in case of the 5-norm to the center of Π.

x x x x f(x) f(x) f(x) f(x) (a) (b) (c) (d)
Remark III.6. Every TDA must show a high sensitivity near the boundary of Π. This is illustrated with the minimal 2norm and the barycenter by the following limit case: with the same Π as before, let W(t) = (1 -t, -1) and w e = 100 + 45|1 -t| -55(1 -t). The time varying linear subspace W(t)τ + w e (t) = 0 is displayed in Figure 3-(a). It intersects only the boundary of Π, hence not satisfying Stater's condition at any time. More precisely, one can observe that for t < 1 the only feasible tension is τ = (10, 100) therefore, both TDAs Figs. 3-(b)(c) will output this tension while, for t > 1, the only feasible tension is τ = (100, 100) and then, both TDAs will output this tension. For t = 1 the whole upper face of Π is feasible so the minimal 2-norm is attained for τ = (10, 110) while the barycenter is τ = (55, 110). The profile τ 1 (t) for both TDAs are depicted in black in Figures 3-(b)(c), and are seen to be discontinuous 5 . Now, enlarging slightly Π as Π = τ | 0 < 9 ≤ τ ≤ 111 , Slater's condition is satisfied at all time and, as expected, the profiles become continuous as shown in gray in the same graphics. The Analytic Center can now be computed and the tension profile of the first cable is depicted in Figure 3-(d). One can see that because of the proximity of the boundary where the tension profile was discontinuous, all TDAs give rise to abrupt change of the tension.

IV. THE ANALYTIC CENTRE TENSION DISTRIBUTION ALGORITHM

A. Definition and main properties

A cable tension set Ξ, defined as the intersection between Π and the set of convex inequality constraints C = {τ | g i (τ ) ≤ 0}, Ξ = Π ∩ C is considered. Then, the (weighted) Analytic Centre is the optimal tension vector τ that minimizes the function ϕ(τ ) defined by

- m i=1 c i log(τ i -τ i ) -c i log(τ i -τ i ) - p i=1 ci log(-g i (τ )), (14)
subject to the equality constraints Wτ + w e = 0. As said previously, the constraints τ ∈ Π are implicit in this problem formulation and are enforced by the objective logarithmic barriers. From now on, we assume that the equality constraints are feasible with respect to these strict inequality constraints, which correspond to the Slater's constraint qualification for convex optimization. In practice, this means that the pose is not on the boundary of the WFW.

Remark IV.1. The requirement of Slater's condition to hold is not a drawback of the Analytic Centre. Most TDAs will fail computing tensions if the Slater condition is not satisfied. Moreover, such poses being exactly on the boundary of the WFW have no practical interest.

The logarithmic barriers are strictly convex constraints in all direction which, together with the convexity of logarithmic barriers of the nonlinear constraints, make the function ϕ(τ ) strictly convex. Therefore the Analytic Centre is uniquely defined.

The behavior of TDA for time varying poses is of critical importance. Discontinuities in tension profiles or in their derivatives can create unwanted vibrations on the end-effector. Consequently, one asks for tension profiles to be as smooth as possible. Formally, suppose that the pose p(t) and the external wrench w e (t) are k times differentiable. The expression of the wrench matrix W(p(t)) shows it is as smooth as p(t). To lighten the notation, in the rest of the paper, writing W(t) is as W(p(t)). Then the Analytic Centre τ (t) is uniquely defined for each time instant by

τ (t) = argmin W(t)τ +we(t)=0 ϕ(τ ). (15)
The following Theorem IV.1 shows the smoothness of τ (t).

A proof relying on the application of the implicit function theorem for k times differentiable function is possible because the optimization problem Eq.( 15) is strictly convex with no inequality constraints 6 .

Theorem IV.1. Let's assume the cost function ϕ(τ ) strictly convex and, at each instant of time, the optimization problem feasible 7 while W(t) is full rank. Provided that W(t) and w e (t) are k 1 ≥ 0 times differentiable and that ϕ(τ ) and g i (τ ) are k 2 ≥ 2 times differentiable, then the unique solution τ (t) of the time dependent optimization problem Eq.( 15) is min{k 1 , k 2 -1} times differentiable.

Proof. Let us consider an arbitrary time instant t 0 and define τ 0 = τ (t 0 ) the unique solution of the optimization problem Eq.( 15) for W 0 = W(t 0 ), w e0 = w e (t 0 ). Since the implicit inequality constraints τ < τ < τ and g i (τ ) < 0 are strict and have, by assumption, some common solutions with the linear constraints W(t)τ + w e (t) = 0, they cannot be active and Lagrange theorem can be applied to show that the solution τ 0 of the optimization problem (15) satisfies the system of equations

∇ϕ(τ ) + W T λ = 0 and W τ + w e = 0 (16)
for W = W 0 and w e = w e0 , where λ is the vector of Lagrange multipliers. This is a square system of equations with variables τ ∈ R m , λ ∈ R n and where W and w e are regarded as parameters. Every variables and parameters appear linearly in the system Eq.( 16). Note that ∇ϕ(τ ) is k 2 -1 ≥ 1 times differentiable with respect to τ (because ϕ(τ ) is supposed k 2 ≥ 2 times differentiable). Therefore, the whole system Eq.( 16) is k 2 -1 times differentiable with respect to all variables and parameters. Its Jacobian matrix with respect to variables τ and λ is

∇ 2 ϕ(τ ) W T W 0 n×n , (17)
6 Non strict convexity or change in constraint activation entails nonregularities of the first order conditions Jacobian matrix that prevent applying the implicit function theorem as in the proof of Theorem IV.1.

7 I.e., the cable tension set defined by τ < τ < τ and g i (τ ) < 0 intersects the linear subspace W(t)τ + we(t) = 0.

which is well known to be regular provided that ∇ 2 ϕ(τ ) is regular and that W is full rank. Those conditions are satisfied for τ 0 (because ∇ 2 ϕ(τ 0 ) is symmetric definite positive due to the strict convexity of ϕ(τ )) and W 0 . Therefore, one can apply the implicit function theorem to characterize the variation of the solution τ 0 when parameters W 0 and w e0 change: there exist neighborhoods N W of W 0 and N we of we 0 and solutions τ * (W, w e ) and λ * (W, w e ) defined inside those neighborhoods and k 2 -1 times differentiable, which satisfy the system Eq.( 16). Thus, since W(t) and w e (t) are continuous, by definition of the continuity there exist a neighborhood N t of t 0 such that, for all t ∈ N t , W(t) ∈ N W and w e (t) ∈ N we . In these neighborhoods, the function τ * (W, w e ) is k 2 -1 times differentiable while functions W(t) and w e (t) are k 1 times differentiable by assumption. Therefore, their composition τ * W(t), w e (t) is min{k 1 , k 2 -1} times differentiable. Finally, since the solution to the system of equations at each time instant is unique, τ (t) = τ * W(t), w e (t) inside N t and it is min{k 1 , k 2 -1} times differentiable at t 0 .

B. Practical computation of the Analytic Centre

This section aims to show how the use of a suitable solver, to solve the Analytic Centre optimisation problem Eq.( 15), can reduce the computation time of the solution w.r.t. build-in solvers demostrating real-time capabilities. In particular, the benchmark is made by considering the Sequential Quadratic Program (SQP) algorithm [START_REF] Nocedal | Numerical Optimization[END_REF] as a representative among the existing one in the MATLAB library. The choice lie on SQP since it results to be the fastest among the available ones in solving Eq. [START_REF] Verhoeven | Tension distribution in tendon-based stewart platforms[END_REF]. On the other hand, the Newton's Algorithm 1 is used [START_REF] Boyd | Convex optimization[END_REF]. It is well-known for its quick convergence and sensibility to the choice of initial iteration. This latter can cause several problems in terms of convergence. Anyway some stratagems that make the algorithm more robust exists. Indeed, it is possible to adapt Newton's algorithm to manage infeasible starting points 8 and update the iterant at each solution. Moreover, since the objective function of Analytic Centre is strictly convex it facilitates the convergence to the solution.

All the simulations are performed in MATLAB using a MacBook Pro Retina 2015 with an Intel Core i5 2.7 GHz processor and 8 GB RAM 1867 MHz DDR3.

With this aim in mind, let us compare the computational cost of Newton and MATLAB fmincon-SQP methods to solve the minimization problem in the form Eq. ( 15). The scope are to extract and compare the iteration number and computational time of mentioned methods while verifying that the obtained tension profiles coincide. Hereby, the study case considered resolves in a CDPR with 4 cables as depicted in Figure 5-(a), the point-mass is supposed to follow a circular trajectory Eq. ( 18) of radius r c = 0.5 m in 10 seconds where velocity and acceleration are null at the start and end of the path. The mass of the load is set to be 20 kg and it is guided by four cables whose tension limits are fixed to τ = 50 N and τ = 400 N , respectively.

x(t) = r c (cos(2πs(t)) + 1.75) s(t) ∈ [0, 1] y(t) = r c (sin(2πs(t)) + 1.75) t ∈ [0, 10]s [START_REF] Hassan | Optimization of actuator forces in cablebased parallel manipulators using convex analysis[END_REF] with s(t) being a 7 th -degree polynomial with abovementioned boundary condition acting up to the 3 rd derivative while the objective function is

ϕ(τ ) = - 4 i=1 log(τ -τ i ) + log(τ i -τ ). ( 19
)
As far as the computational performances are concerned, the data acquired during the simulation, useful for the comparison, are collected in Table I. Moreover, in Fig. 4, the computational time necessary to find a solution for each time step' simulation is reported. Note that at each call of the solver, the initial iterant is updated with the previous solution found. Further, observe that the gradient of the cost function ϕ(τ ) was also provided. Both these expedients are employed to speed up the computation.

Therefore, analysing the data, it seems clear that implementing an ad-hoc algorithm to solve this optimization problem Eq. ( 15) helps reducing the computational cost w.r.t. buildin functions. Indeed, Newton is generally 2-order faster than SQP.

Note coherency between the maximum number of iteration and first call of the Newton method. With SQP, the maximum iteration number does not coincide with the first call. This discrepancies is due to the different architecture of the two algorithms which tackle the problem differently. Briefly, the SQP tries to approximate the objective with a quadratic model and solve the sequence of subproblems. Hence, it generally requires more iterations (and time) to converge. Though the maximum and the initial number of iteration does not coincide for the SQP, the computational time required for the first iteration result coherent, for both algorithms, as the most expensive due to the distance between the initial (i.e. tentative) and the first solution.

To conclude, is worth noticing the Newton trend in Fig. 4: it is evident the presence of two steps at the start and end of the simulation. These are due to non-homogeneous distribution of the points along the trajectory which, because of s(t), are more dense about the initial and final of the path. The quick convergence of the method makes it sensible to the discretization of the trajectory.

Moreover, looking at both trends, it seems that Newton oscillates more than SQP. This effects is due to passive processes ongoing on the laptop. Indeed, their influence affect more the Newton trend since it works 2-order faster than SQP: a disturbance has minor effect over a big quantity.

V. SIMULATIONS AND CASE STUDIES

The implementation of the Analytic Centre, defined in Sec. IV-C, to different case studies are conveyed in this section. In accordance with previous tension distribution methods, Step time comparisons between Newton and fmincon-SQP algorithms. This graph represents the mean value of computational times. It is averaged over three simulations in order to reduce the passive effects of other processes running in parallel on the processor.

mostly developed for CDPRs, this section will follow the same explanatory line to emphasize the peculiarities of the Analytic Centre, even though it applies for a wide range of CDPs [START_REF] Jamshidifar | Static workspace optimization of aerial cable towed robots with land-fixed winches[END_REF], [START_REF] Masone | Cooperative transportation of a payload using quadrotors: A reconfigurable cable-driven parallel robot[END_REF]. In particular, planar CDPRs with four and six cables and a planar MCDPR with four cables will be considered.

In addition, also a spatial CDPR with eight cables is studied. Their dynamic will be simulated while cable tension profiles are computed using the Analytic Centre and various stateof-the-art methods. Indeed, the purpose of this section is to show and compare the capabilities of the Analytic Centre with respect to previous works in order to get its pros and cons. With this in mind, a total of four examples will be considered to reveal its main peculiarities.

In the first case study, Sec. V-A, the comparisons with existing methods will focus on the robustness, continuity and differentiability of the solutions.

The second example, Sec. V-B, aims to highlight the generality of the proposed method illustrating the possibility to apply it with any DoR while providing feasible, robust and smooth tension profile solutions. Furthermore, the computa- tional time of the methods that are capable of providing a feasible and continuous solution with general DoR is compared to confirm that the proposed method can operate in real-time.

Algorithm 1 Newton Algorithm Input: τ 0 % generally τ While: ||Eq.(16)|| > Tol & iter < Max iter Compute: τ k+1 Update: τ 0 → τ k+1 End Output: τ % last τ 0 A A 1 2 A3 A4 A1 A2 A3 A4 A5 A6
The third case study described in Sec. V-C lends itself well to emphasize the versatility of the Analytic Centre. Indeed, it shows that it is possible to take into account non-linear timevarying constraints in a natural way returning, again, a set of tensions with the aforementioned characteristics.

To conclude, the last simulation, Sec. V-D, somehow summarizes the previous results Analytic Centre capabilities even when a spatial case is considered.

A. Planar CDPR composed of four cables and a point-mass end-effector

The architecture of the robot considered in this study as well as the trajectory is the same used in Section IV-B.

During the tracking task, the robust index, defined in Eq.( 10) is computed. Moreover, several state-of-the-art methods are reproduced to determine the tension profiles according to different criteria. The obtained results are collected and depicted in Figure 6. At first sight, all implemented methods provide feasible solutions: the tension profiles found satisfy the maximum and minimum tension constraints Eq.( 4) (represented in each figure by the dashed red line) and the load equilibrium Eq.( 1). Moreover, the profiles look not symmetric since the trajectory has been shifted in the left-bottom corner of the robot workspace. Reproduced methods are now investigated one-byone in order to remark on the relevant aspects necessary for the comparison.

Let's consider the QP method first whose tension profiles are depicted in Fig. 6-(d) [START_REF] Taghirad | An analytic-iterative redundancy resolution scheme for cable-driven redundant parallel manipulators[END_REF]. This approach tends to save energy required from the actuators by minimizing the tension components as much as possible. Indeed, the results show that two tension components τ 1 and τ 2 (cables below the load) are found to assume the minimum value. As a result, the robustness index is the lowest among the other methods as shown in Fig. 6-(f). In practice, tension profiles are continuous and differentiable but the robot suffers from low stiffness which can lead to undesired vibrations and oscillations of the platform.

On the other hand, the Robust solution reported in Fig. 6-(c) is the one that has the highest robust index value. Anyway, since this approach belongs to the family of LP optimization methods [START_REF] Borgstrom | Rapid computation of optimally safe tension distributions for parallel cable-driven robots[END_REF], the solution can result in discontinuities. This example reveals discontinuities of the solution while tracking a smooth circular trajectory. In some sense, achieving the maximum stiffness is akin admitting discontinuities in tension profiles. Therefore, the presence of discontinuities, even when continuous trajectory and wrench are considered, can cause unwanted phenomena (e.g. vibrations) on the end-effector.

In Figure 6-(b) the solution found using the Barycentric approach [START_REF] Mikelsons | A real-time capable force calculation algorithm for redundant tendon-based parallel manipulators[END_REF], [START_REF] Gouttefarde | A versatile tension distribution algorithm for n-dof parallel robots driven by n + 2 cables[END_REF] is similar in both tension profiles and robustness to the Analytic Centre solution reported in Fig. 6-(a). Despite their affinity, there are crucial differences among them. Indeed, although the Barycentric provides robust and continuous profiles there is no formal evidence of their differentiability. Moreover, this technique is not general since is valid for DoR ≤ 2. This latter will be better clarified in the next section.

The Improved-Closed Form [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF] approach whose solution is displayed in Fig. 6-(e) furnishes continuous tension profiles with a close, but generally lower, robust index value than the Analytic Centre method. It constitutes an improvement to the previous work [START_REF] Pott | Closed-form force distribution for parallel wire robots[END_REF]. Anyway, the absence of a formal proof guaranteeing the convergence to a feasible solution represents the main shortcoming of this approach.

So far, for mentioned literature methods the robustness, continuity and differentiability of the tensions were discussed. To sum up, these examples showed that, among these TDAs, the Analytic Centre is the only one able to preserve mentioned properties simultaneously.

B. Planar CDPR composed of six cables and a rigid-body platform end-effector

The use of several cables serves as an example to demonstrate the generality w.r.t. DoR of the proposed method in finding a robust solution which is also continuous and differentiable. The scheme of the robot is depicted in figure 5-(b), the centre of mass of the square platform is supposed to follow the same trajectory as in the previous case Sec.V-A. The axes of the mobile reference frame will remain parallel to those of the fixed ones throughout the simulation. The mass of the load is 65 kg and it is guided by six cables whose tension limits are the same as in the previous example. The choice of a large mass forces the system to the boundary of the WFW.

Simulation results are given in Figure 7. Observe that, because of several cables DoR = 3 and therefore the Barycentric approach can not be applied. Indeed, the existing method [START_REF] Mikelsons | A real-time capable force calculation algorithm for redundant tendon-based parallel manipulators[END_REF], [START_REF] Gouttefarde | A versatile tension distribution algorithm for n-dof parallel robots driven by n + 2 cables[END_REF] relies on finding the barycenter of a 2D polytope Λ, defined in Eq. [START_REF] Surdilovic | String-man: Wire-robot technology for safe, flexible and human-friendly gait rehabilitation[END_REF], by means of triangulation techniques. The latter can be efficiently applied in low-dimensional space. Already from a dimension higher than two this technique results impractical due to the increase in its computational cost. Mentioned issue combined with the common structure of CDPRs, which often foresees a maximum of DoR = 2, prevented researchers from a generalization. Anyway, the advent of CDP where a large redundancy can occur is, among the others, another motivation for the introduction of the Analytic Centre.

The Improved Closed-Form method instead fails furnishing a feasible solution. This happens because the high mass pushes the cable tensions outside the limits. Hence, although it is an improved version of the Closed-Form method [START_REF] Pott | Closed-form force distribution for parallel wire robots[END_REF], it is one is not always able to converge in DoR steps. The same authors express this possibility in their work [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF]. Thus, this case study, identifies a reduced WFW coverage of the Impoved Closed-Form w.r.t. the Analytic Centre, QP and Robust methods which find a feasible solution.

Among the remaining methods, whose solutions are reported in Fig. 7, the Analytic Centre approach appears to be the best compromise between robustness and continuity of the solution. Indeed, looking at Figure 7-(d), the QP maintains the lowest robust index reducing the energy consumption at the price of a reduced stiffness of the robot whereas the Robust method registers the highest index value but introduces discontinuities of tensions in time which can cause unwanted phenomena as well.

However, at this point, someone might wonder what are the performance of this method in terms of computational time. Hence, to make a fair comparison, the computational time required to determine the solution to each optimization problems is considered. In this regard, a comparison of the approaches that return regular continuous and feasible tension profiles for any DoR is given.

The graph in Fig. 8 shows the time required by the solver to find a solution during the simulation for both Analytic Centre and QP methods. In particular, these data are generated by averaging the solution times over 3 simulations to remove disturbances that can occur due to passive processes ongoing on the computer. For the Analytic Centre approach the Newton algorithm discussed in IV-B is used to solve the optimization problem whereas for QP method, the Sequential Quadratic Programming (SQP) algorithm [START_REF] Nocedal | Numerical Optimization[END_REF] is a natural choice given the nature of the problem. The gradient of the objective functions is provided in advance as well as the iterant is updated at each call in order to speed up the solvers.

Looking at the graph Fig. 8, the computational times of the Newton methods applied to the Analytic Centre is 2order smaller than the QP. The differences in the iteration time are to be attributed to the diversity between the solvers, the objective functions and the absence of active constraints (i.e. inequalities) in the case of the Analytic Centre approach. The inequalities which define the QP optimization problem increase its complexity forcing the solver consuming time for the few iterations required to converge. In both cases, the first iteration is the one which costs more. The maximum number of iterations recorded amounts to n iter,AC = 7 and n iter,QP = 3. Because the QP method is generally applied in real-time application, it is reasonable to argue that also the Analytic Centre can be applied in real-time to control cable the tensions of a CDP. Table II summarizes the discussed data used for time comparisons.

C. Planar MCDPR composed of four cables and a point-mass end-effector

This example takes its cue from mobile cable systems [START_REF] Tahir Rasheed | Tension distribution algorithm for planar mobile cable-driven parallel robots[END_REF], [START_REF] Rasheed | Wrench-feasible workspace of mobile cable-driven parallel robots[END_REF]. Indeed, in these systems, the choice of the cable tensions has to take into account the static friction between the ground and the Mobile Base (MB) to avoid sliding effects during the execution of a task. This motivates the need to introduce non-linear constraints in the cable tensions computation. The scheme of the MCDPR considered in this example is given in Fig. 9. It constitutes a simplified version of the original MCDPRs. Indeed, in practice, friction is considered to act on the four wheels that move each MB instead of assuming its influence as concentrated on its centre of mass. However, this architecture is enough for the scope of this section and although it constitutes a simplified model, it does not represent a limitation for the presented method. Later on, the MBs will be considered as fixed at the ground and therefore the robot will coincide with the study case discussed in Sec. V-A. This latter will allow to make comparisons and see the consequences on the MCDPR due to the differences in the tensions distribution when friction is not considered.

For the current case, the mass of the moving bases is considered to be equal to m MB1 = m MB2 = 65kg while the load mass is 5kg. The static friction coefficient is µ = 0.5.

A A

1 2

A3 A4

Centre of mass The load trajectory is the same as in the previous cases as the tension limits. The equations which constitute the non-linear constraints C for the cable tensions are the followings

||τ 1 u 1 + τ 4 u 4 || 2 2 ≤ m 2 MB1 g 2 µ 2 , (20a) ||τ 2 u 2 + τ 3 u 3 || 2 2 ≤ m 2 MB2 g 2 µ 2 , (20b)
these can be derived by using the free-body diagram for both MB 1 and MB 2 . Indeed, the MB can be considered as a decoupled systems.

Results of the computed tension profiles are depicted in Fig. 10. The objective function used for this simulation is reported in Eq. ( 21) for sake of clarity.

ϕ(τ ) = - 4 i=1 log(τ i -τ min,i ) + 4 i=1 log(τ max,i -τ i )+ log -τ 1 + τ 4 2 + m 2 MB1 g 2 µ 2 + log -τ 2 + τ 3 ) 2 + m 2 MB2 2 µ 2 . ( 21
)
The tension profiles show how the introduction of friction narrows the set tensions that can be used to perform a task. In particular, the tensions, Fig. 10 continuous line, assume reduced values compared to the case of fixed moving bases. In other words, avoiding sliding conditions for the MBs means reducing tension values while performing a task. The consequences of neglecting friction while computing tension profiles can be visualized in Fig. 11. The trend of the constraints demonstrates that neglecting friction can cause the MBs to slide since the dotted curves exceed the static friction limit. An alternative way to display the same phenomena consists in investigating if the λ vectors of the two solutions, for t = 3s and t = 5s, belong to the intersection between the feasible polygon and the non-linear constraints Ξ = Λ ∩ C; the dual of Ξ. The equations necessary to map the constraints into Λ, as done in Figs. 12 and 13, are the followings Fig. 10. Cable tension results for the MCDPR: continuous profiles refer to the case with non-linear constraints whereas dotted ones do not consider any additional constraints (i.e. the MBs are considered as fixed). 20) for the two moving body. Constraints gτ 1 τ 4 and gτ 2 τ 3 take friction into account while gτ 1 τ 4 and gτ 2 τ 3 do not. It can be seen that for the present case, the gτ i τ j (dotted profiles) constraints exceed the sliding limit several times.

τ p,1 + (Nλ) 2 1 + 2τ p,1 (Nλ) 1 + τ 2 p,4 + (Nλ) 2 4 + 2τ p,4 (Nλ) 4 + 2 τ p,1 τ p, 4 + τ 1,p (Nλ) 4 + τ 4,p (Nλ) 1 + (Nλ) 1 (Nλ) 4 ≤ m 2 MB1 g 2 µ 2 , (22a) τ 2 p,2 + (Nλ) 2 2 + 2τ p,2 (Nλ) 2 + τ 2 p,3 + (Nλ) 2 3 + 2τ p,3 (Nλ) 3 + 2 τ p,2 τ p, 3 + τ 2,p (Nλ) 3 + τ 3,p (Nλ) 2 + (Nλ) 2 (Nλ) 3 ≤ m 2 MB1 g 2 µ 2 , (22b)
where τ p,i and (Nλ) i represent the components of the particular and general solution presented in Eq. ( 6), respectively. Thus, the constraint Equations ( 22) in Λ can be retrieved by substituting Eq. ( 6) into Eqs. [START_REF] Mikelsons | A real-time capable force calculation algorithm for redundant tendon-based parallel manipulators[END_REF]. In other words, it suffices to employ the composition rule for functions. Further, observe that the elements (Nλ) i could be thought as τ g,i but, here, the explicit version with λ i is clearer since the aim is to investigate the conic' constraints intersection with the feasible polygon Λ. Indeed, as a confirm of what displayed in Fig. 11, Figs. 12 and 13 demonstrate that the sliding condition occurs when λ / ∈ Ξ. This example showed the ability of the Analytic Centre to deal with non-linear constraints, so far approximated as linear as in [START_REF] Tahir Rasheed | Tension distribution algorithm for planar mobile cable-driven parallel robots[END_REF], [START_REF] Rasheed | Wrench-feasible workspace of mobile cable-driven parallel robots[END_REF], for example. In addition, as done in Eq. ( 21), the possibility to superimpose general (linear and non-linear) time-varying constraints on specific cable tensions makes this approach compatible with Human Physical-Interaction (HP-I) or collaborative applications. Indeed, mentioned applications A1 A2 A3 A4 A5 A6 A7 A8 Fig. 14. Spatial architecture considered, see [START_REF] Métillon | Performance and interaction quality variations of a collaborative cable-driven parallel robot[END_REF] for details. ask to cope with humans thus demanding taking precautions in terms of tension limit along specific cables which, usually, are close to the human [START_REF] Métillon | Performance and interaction quality variations of a collaborative cable-driven parallel robot[END_REF].

D. Spatial CDPR composed of eight cables and a rigid-body end-effector

This example is intended to summarise the distinctive features of the TDAs analysed and the one proposed here. In particular, a spatial CDPR designed for collaborative tasks will be considered [START_REF] Métillon | Performance and interaction quality variations of a collaborative cable-driven parallel robot[END_REF]. This robot considers eight cables and has 2 DoRs. The cables are fixed at the upper part of the frame, near the vertices of the cubic cell, as shown in the Figure 14 9 . This choice avoids possible dangerous contact with the operator intent on collaborating with the platform.

For the purposes of the paper, it is sufficient to consider a tracking task: the trajectory under consideration is always circular and arranged in the horizontal plane Eq.( 23), similar to the previous cases. The full dimension of the robot and its precise description can be found in [START_REF] Métillon | Performance and interaction quality variations of a collaborative cable-driven parallel robot[END_REF]. The mass of the platform is m = 7kg and it is guided by eight cables whose tension limits are fixed to τ = 1 N and τ = 100 N , respectively.

    

x(t) = r c cos(2πs(t)) s(t) ∈ [0, 1] y(t) = r c sin(2πs(t)) t ∈ [0, 10]s z(t) = 2 [START_REF] Bouchard | On the ability of a cabledriven robot to generate a prescribed set of wrenches[END_REF] As done for the other case study, several TDAs are compared. Their tension profiles and robustness indices are reported in Fig. 15. Alike the case study in Section V-A all the methods provide feasible tension profiles and pertain to the main characteristics discussed above. However, as seen before, existing methods results in optimal performances only for particular predefined operating conditions. In fact, in many cases (for certain tasks), some of their relevant properties are lost such as of tension profile, quick convergence, generality for any DoR and so on. 9 The real prototype is located in Nantes, France at LS2N.

A further aspect to be examined and added to the previous considerations is highlighted by the architecture considered. Indeed, this CDPR is redundant but not fully-constrained as the load is suspended. In the several case studies analysed in the state of the art, this type of architecture is not very common but nevertheless important for recent applications.

In particular, the major difference that is evident, compared to the previous examples, concerns the values taken by the robustness indices.

In fact, this time, the QP method, see Fig. 15-(d), does not have the lowest index. Moreover, the other methods show relatively little difference in terms of robustness. In some ways, they can be considered equivalent.

When it is the case of suspended loads, it has been seen [START_REF] Paola | A preliminary study of factors influencing the stiffness of aerial cable towed systems[END_REF] that the capability to optimise the cable tension is reduced. In other words, the absence of cables under the platform leads to a reduction of the area (hyper volume, in general) of the polygon Λ and thus in a reduced ability to generate a wrench on the load. This also explains why the range of variation of the various robustness indices is small compared to previous cases (all fully-constrained). All indices are very similar to each other with the robust method, as always, that holds the maximum value. To conclude, this example confirm that, in general, existing methods results in optimal performances only for particular predefined operating conditions (i.e. task at hand and architecture). Instead, the Analytical Centre method always guarantee the best compromise in terms of robustness, smoothness of tension profiles, generality and reliability.

VI. CONCLUSIONS

In this paper, the tension distribution problem for Cable-Driven Platforms (CDPs) with DoR ≥ 1 was addressed. The formulation of the optimization problem with the barrier function enabled eliminating inequality constraints, considering non-linear constraints and reaching a robust and unique solution in tension space. Therefore, this made unnecessary to develop additional algorithms for building Λ polytope that are complex for high Degree of Redundancy (DoR). Moreover, using the Analytic Centre, cable tension profiles were proven to be continuous and differentiable. Proposed examples showed that the Analytic Centre is the best compromise among reproduced methods in terms of robustness, continuity and differentiability. The computational aspect was also considered, in fact, the implementation of the Newton method allowed to demonstrate its real-time capabilities. Finally, this method is not only general, but its versatility extends its applicability to a wide range of cable robots such as MCDPRs and ACTSs. Therefore, it lays the first brick towards tension management for collaborative tasks. Future work deals with the implementation of this technique on a real framework.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Barrier function: (a) the cost function f (x) = x subject to g(x) = x(x -1) ≤ 0. Remaining graphics (b),(c) and (d): the unconstrained problem fc(x) = c log(-x(x -1)) with c ∈ {1, 0.1, 0.01}. One can observe that the unconstrained minimizer converges toward the constrained minimizer as c converges to 0.

Fig. 4 .

 4 Fig. 4. Step time comparisons between Newton and fmincon-SQP algorithms. This graph represents the mean value of computational times. It is averaged over three simulations in order to reduce the passive effects of other processes running in parallel on the processor.

Fig. 5 .

 5 Fig. 5. Two planar CDPRs: (a) composed of four cables and a point-mass end-effector (b) composed of six cables and a rigid-body platform end-effector

Fig. 6 .

 6 Fig. 6. Tension profiles for: (a) Analytic Centre approach (b) Barycentric method (c) Robust (LP) technique (d) Quadratic Programming (QP) (e) Improved-Closed Form. The trends of the robust index for mentioned methods are reported in Figure (f).

Fig. 7 .

 7 Fig. 7. Tension profiles for: (a) Analytic Centre (b) Robust (LP) method (c) Quadratic Programming (QP). The trends of the robust index for mentioned methods are reported in Figure (d).

Fig. 9 .

 9 Fig. 9. Architecture of a planar MCDPR with 4 cables

2 Fig. 11 .

 211 Fig. 11. Trend of non-linear constraint values Eq.(20) for the two moving body. Constraints gτ 1 τ 4 and gτ 2 τ 3 take friction into account while gτ 1 τ 4 and gτ 2 τ 3 do not. It can be seen that for the present case, the gτ i τ j (dotted profiles) constraints exceed the sliding limit several times.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Intersection between feasible polygon Λ and the non-linear constraints Eq. 20. This plot depicts the case at t = 3s where both λ vectors are feasible.

Fig. 15 .

 15 Fig. 15. Tension profiles for: (a) Analytic Centre approach (b) Barycentric method (c) Robust (LP) technique (d) Quadratic Programming (QP) (e) Improved-Closed Form. The trends of the robust index for mentioned methods are reported in Figure (f). Observe that the colors used for the tension profiles coincide with the colors of the cables in Fig. 14.

  Computational time per solution during the simulation. This graph represents the mean value of computational times. It is averaged over three simulations in order to reduce the passive effects of other processes running in parallel on the processor.

		TABLE II
	SUMMARY OF THE SIMULATION DATA
	Results	Newton (AC)	fmincon-SQP (QP)
	Time Step	0.001s	0.001s
	n iter first Solution Mean Computational Time Max. Iteration Number	7 7.3e -5 s sol. 7	2 6.8e -3 s sol. 3
	Stopping Criterion	||Eq.(16)|| 2 < 10 -10	fmincon default options
	Fig. 8.

DoR = m -n where m is the cable number and n stands for the Degree of Freedom (DoF ) of the platform.

Both TDAs are not well defined in case of cable dependent tensions limits since the aim is to maximize the distance d(τ ) from the boundaries τ and τ .

Or equivalently, the distance to the closest face of Π.

This condition is typical in convex optimization and is called Salter's condition. Data not satisfying this condition have no Analytic Centre. This correspond to the very boundary of the WFW, and has no practical impact.

The statement and the proof of the continuity of the Barycenter presented in[START_REF] Mikelsons | A real-time capable force calculation algorithm for redundant tendon-based parallel manipulators[END_REF] do not involve Slater's condition, so the example provided here shows it must be not fully correct. The proof of continuity of the norm based TDA given in[START_REF] Verhoeven | Analysis of the workspace of tendon-based stewart platforms[END_REF] is difficult, involving continuity of set-valued mappings and Berge maximum theorem. Condition (4.53) in Lemma 4.6 requires that the linear subspace has full dimension within Π, which seems closely related to Slater's condition.

This means that the initial iterant satisfies the tension limits but does not necessarily satisfy the equilibrium Eq.(1).

Vincenzo Di Paola received the B.Sc. and M.Sc. in mechanical engineering from the University of Bologna with a focus on automation and robotics in 2017 and 2020, respectively. He is currently pursuing a double-degree Ph.D. in robotics and mechatronics with the University of Genova, Italy and the Ecole Centrale de Nantes, France. After studying persistent screw systems during his master's thesis, in 2020, he worked as a Research Fellow for three months with the Department of Industrial Engineering (DIN), University of Bologna, where he started his actual research on the design and control of cable-driven parallel robots, aerial cable-towed systems and their hybrid versions.


Lucie Leboulleux Abc , Mamadou N'diaye

Eldorado Riggs

Egron Sylvain

Johan Abc

Mazoyer

Laurent Pueyo

Elodie Choquet

Marshall D Perrin

Jeremy Kasdin

Jean-François Sauvage

Thierry Fusco

Rémi Soummer

High-contrast imager for Complex Aperture Telescopes (HiCAT). 4. Status and wavefront control development

Keywords: exoplanets, high-contrast imaging, wavefront sensing, wavefront control, Speckle Nulling, vibration analysis

Segmented telescopes are a possible approach to enable large-aperture space telescopes for the direct imaging and spectroscopy of habitable worlds. However, the increased complexity of their aperture geometry, due to their central obstruction, support structures and segment gaps, makes high-contrast imaging very challenging. The High-contrast imager for Complex Aperture Telescopes (HiCAT) was designed to study and develop solutions for such telescope pupils using wavefront control and starlight suppression. The testbed design has the flexibility to enable studies with increasing complexity for telescope aperture geometries starting with off-axis telescopes, then on-axis telescopes with central obstruction and support structures (e.g. the Wide Field Infrared Survey Telescope [WFIRST]), up to on-axis segmented telescopes e.g. including various concepts for a Large UV, Optical, IR telescope (LUVOIR), such as the High Definition Space Telescope (HDST). We completed optical alignment in the summer of 2014 and a first deformable mirror was successfully integrated in the testbed, with a total wavefront error of 13nm RMS over a 18mm diameter circular pupil in open loop. HiCAT will also be provided with a segmented mirror conjugated with a shaped pupil representing the HDST configuration, to directly study wavefront control in the presence of segment gaps, central obstruction and spider. We recently applied a focal plane wavefront control method combined with a classical Lyot coronagraph on HiCAT, and we found limitations on contrast performance due to vibration effect. In this communication, we analyze this instability and study its impact on the performance of wavefront control algorithms. We present our Speckle Nulling code to control and correct for wavefront errors both in simulation mode and on testbed mode. This routine is first tested in simulation mode without instability to validate our code. We then add simulated vibrations to study the degradation of contrast performance in the presence of these effects.

INTRODUCTION

The next generation of space telescopes for direct imaging and spectroscopy of exoplanets includes telescopes with a monolithic mirror, such as the Wide Field Infrared Survey Telescope (WFIRST), [START_REF] Spergel | Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report[END_REF] and Large Ultra-Violet Optical Infrared (LUVOIR) telescopes with segmented primary mirror, like ATLAST [START_REF] Postman | Advanced Technology Large-Aperture Space Telescope: science drivers and technology developments[END_REF][START_REF] Feinberg | A Cost-effective and Serviceable ATLAST 9.2m Telescope Architecture[END_REF] or HDST. [START_REF] Dalcanton | From Cosmic Birth to Living Earths: The Future of UVOIR Space Astronomy[END_REF] Because of the complexity of their pupils, high-contrast imaging becomes more challenging, including starlight suppression with coronagraphy and wavefront sensing and wavefront control for contrast stability in the presence of vibrations.

The High-contrast imager for Complex Aperture Telescopes (HiCAT) testbed has been developed to enable studies on starlight suppression, wavefront sensing (WFS), and wavefront control (WFC) for such unfriendly pupils. New coronagraph designs for any kind of aperture that are based on the Apodized Pupil Lyot Coronagraph (APLC) [START_REF] N'diaye | Apodized Pupil Lyot Coronagraphs for Arbitrary Apertures. IV. Reduced Inner Working Angle and Increased Robustness to Low-order Aberrations[END_REF][START_REF] N'diaye | Apodized Pupil Lyot Coronagraphs for Arbitrary Apertures. V. Hybrid Shaped Pupil Designs for Imaging Earth-like planets with Future Space Observatories[END_REF] are currently developed on simulation and will then be implemented on the testbed. Furthermore, new high-contrast methods will be implemented and tested, such as pupil remapping techniques (e.g. Active Control of Aperture Discontinuities [ACAD] [START_REF] Pueyo | High-contrast Imaging with an Arbitrary Aperture: Active Compensation of Aperture Discontinuities[END_REF][START_REF] Mazoyer | Correcting for the effects of pupil discontinuities with the ACAD method[END_REF][START_REF] Mazoyer | Active compensation of aperture discontinuities for WFIRST-AFTA: analytical and numerical comparison of propagation methods and preliminary results with a WFIRST-AFTA-like pupil[END_REF] ) that use two deformable mirrors (DMs) to convert complex pupils into friendly apertures for coronagraphy.

In previous papers, we presented the overall goals and design of HiCAT [START_REF] N'diaye | High-contrast imager for Complex Aperture Telescopes (HiCAT): 1. Testbed design[END_REF][START_REF] N'diaye | High-contrast Imager for Complex Aperture Telescopes (HICAT): II. Design overview and first light results[END_REF] and described its assembly and integration including the installation and first tests with single Boston Micromachines deformable mirror (DM). [START_REF] N'diaye | High-contrast imager for complex aperture telescopes (HiCAT): 3. first lab results with wavefront control[END_REF] In this communication, we report on progress and challenges since then. In particular we describe the detection and analysis of instability that prevents us from implementing high-contrast imaging codes on HiCAT. From these studies, we derive a new DM mount design and the implementation of a Speckle Nulling code in simulation, without and with instability.

In Section 2, we present the main features of our testbed and our first wavefront control tests. We show that the preliminary results did not reach the desired contrast levels because of the instability issues that were detected during our experiments. In Section 3, we analyze the point spread function (PSF) motion that was observed on the testbed and we derive conclusions on the instability issues. Finally, in Section 4, we present the results of our simulations on wavefront control using a single DM and a classical Lyot coronagraph, including instability effects.

HICAT TESTBED IMPLEMENTATION AND STATUS

Optical and opto-mechanical designs

The detailed optical design of the testbed was presented in N'Diaye et al. [START_REF] N'diaye | High-contrast imager for Complex Aperture Telescopes (HiCAT): 1. Testbed design[END_REF] and it will at the end simulate a complex aperture telescope on which the starlight will be suppressed thanks to an APLC and wavefront control, using two DMs. The design goal was to minimize amplitude-induced errors due to Talbot effects. In addition, the access to several pupil and focal planes enables the correction of phase errors induced by aperture diffraction. As shown in Fig. 1, the HiCAT testbed is provided with a starlight suppression system, based on the APLC architecture. This includes a reflective apodizer, [START_REF] Balasubramanian | High contrast internal and external coronagraph masks produced by various techniques[END_REF] a reflective Focal Plane Mask (FPM), which currently has a 334 micrometer size hole at its center, and a Lyot Stop. Both the Lyot Stop and the FPM are motorized and can be controlled from a computer. The final design of the apodizer is still under development as part of our ongoing coronagraph optimization studies, 6 so it is currently replaced with a high quality flat mirror and the full starlight suppression system is currently equivalent to a Lyot coronagraph. Two DMs (Boston kilo-DM) enable wavefront sensing and control studies on HiCAT. The first one has been set up in a pupil plane in 2014-2015 and allows first studies to be driven. [START_REF] N'diaye | High-contrast imager for complex aperture telescopes (HiCAT): 3. first lab results with wavefront control[END_REF] The second one, out of pupil plane, is not installed yet and is currently replaced by a high quality flat mirror. Furthermore, since HiCAT should be able to simulate a telescope with a complex geometry pupil, it can be featured with a pupil mask similar to the WFIRST pupil, including the spider and the central obstruction. It will also be provided with a segmented mirror (SM), conjugated with the previous complex pupil to simulate an ATLAST/HDST-like primary mirror. This mirror is an Iris-AO SM with 36 segments. At the very end of the testbed, a beam splitter sends the beam to two cameras, CamF and CamP, respectively installed at a focal plane and at a pupil plane.

The testbed is fully aligned with an excellent wavefront quality of 13nm RMS over a 18mm diameter circular pupil [START_REF] N'diaye | High-contrast Imager for Complex Aperture Telescopes (HICAT): II. Design overview and first light results[END_REF] and now includes all the convergent and flat mirrors, the beam launcher, the pupil mask, the first DM (DM1) which is set in a pupil plane, the FPM, the Lyot Stop, and the two cameras. To reduce the stray light, it also has an enclosure, several black screens and laser beam dumps. HiCAT testbed design done with the software Solidworks, the beam is exported from Zemax. The telescope is simulated by a pupil mask, the segmented mirror and off-axis parabolas. The segmented mirror is conjugated with the pupil mask to form a segmented pupil with central obstruction and spider struts. The off-axis parabolas set the telescope aperture. The wavefront control is done with two deformable mirrors. The coronagraph is composed of an apodizer, a focal plane mask and a Lyot stop.

It will also have to be upgraded with the WFIRST-like pupil mask, the SM conjugated with the previous complex pupil to simulate a ATLAST/HDST-like primary mirror, the apodizer that is currently being designed and a second Boston kilo-DM, out of pupil plane.

Wavefront sensing and wavefront control implementation

In 2015, a first Speckle Nulling code was implemented on the HiCAT testbed, using DM1. It enabled the suppression of a speckle in our final image. For more details about this first implementation, please see N'Diaye et al. 2015. [START_REF] N'diaye | High-contrast imager for complex aperture telescopes (HiCAT): 3. first lab results with wavefront control[END_REF] After this first step, we tested a Focal Plane Wavefront Control (FPWC) code 14 but we did not manage to produce a dark hole with the expected contrast because of some Point Spread Function (PSF) drifts between the different exposures.

INSTABILITY ANALYSIS

In this section, we present our instability analysis and introduce a new DM mount design, that should be an important improvement compared to the current one.

PSF motion analysis

The idea here is to characterize the motion effect on recorded videos and, depending on the kind of motion or shape changing we observe, to find both the source of the vibration and its support on the testbed. For example, the motor driver could be the origin of the vibration, but it does not have an real direct influence on the beam, unless an optical component is loose. In the following sections, we present our main tests and the derived conclusions.

The criteria to characterize the motion are: -the amplitude of the motion around an average position as a function of the time -the direction of the motion (the motion could follow only the Y-axis, for example) -the spectrum of the motion. Our tests were done using a Thorlabs CMOS camera, that can take monochrome videos at a frame rate from 25 fps to 250 fps with a 1280 x 1024 pixel detector, where each pixel is a 5.20 µm square. In addition, the small size of the camera (48.6 mm × 44 mm × 25.7 mm) facilitates access to the different planes of study.

HiCAT at rest

As a reference comparison for the following tests, if HiCAT is at rest, so in normal running conditions (hardwares are on and the enclosure is closed), we obtain the results shown in Fig. 2. No main direction of motion could be obtained from the video analysis and only frequencies below 2 Hz seem to appear. Indeed, above 2 Hz, we cannot see any specific peak in the spectrum. Furthermore, we realigned all the frames from a video, where the coronagraph is set, and we noticed after realignment that the PSF was not moving anymore, except that the shape of the core of the coronagraphic PSF was modified with time. This observation made us conclude that there is at least one optical component moving before the FPM.

Which optical components are most unstable?

To find out which optical component could be responsible for the beam direction to change in the testbed, two sets of tests were done.

First, we set the camera at the end of the testbed, at the focal plane and gently tap every optical component that appears on Fig. 1 except for DM1. The idea here is to check if tapping has an influence on the PSF motion. The three main results from these first tests are given in Fig. 3 and tapping every other optical component has no effect on the PSF instability. For these three optics, so DM2, the apodizer and the Iris AO (actually flat mirrors), we do observe an important change of the amplitude and spectrum of the PSF motion.

As a second set of test, we put successively the camera at each of the five focal planes of the testbed and compare the amplitude of the PSF motion normalized by the value of the focal length until this focal plane, which is equivalent of getting a value of the amplitude of the motion angle of the beam instead of the metric amplitude of the PSF. In theses conditions, the amplitude values can be compared even if the camera is set at different focal planes. Without a surprise, the further the camera is from the entrance of the testbed, the worse the amplitude is. We mainly observe the most significant difference between the two focal planes situated right before and right after the two DMs and the mirror O8: the average motion values are multiplied by 3.20 between these two successive focal planes. These three optical components, DM1 and DM2 and O8, then become principal suspects.

As a final test, since our two previous tests were never done on DM1 that is too easily damageable, we built a spare mount, similar to DM1 mount, and set it up in the beam path of the testbed to study its stability as well. We can conclude from the results we obtain that tapping this mount has a huge influence on the PSF motion amplitude (the average value of the motion amplitude is multiplied by 4.87), which means that DM1 is unstable too and can also be considered as a suspect of the PSF motion.

Which source(s) drive(s) the instability?

This test consists in turning off one by one every possible source of motion in the testbed environment, such as the camera fans, the drivers or the air supply and taking videos at the end of the testbed. After analysing the different recorded videos, we reached the conclusion that the motor driver and mainly the camera fans were responsible for the PSF motion. For the cameras turned off, the average motion amplitude indeed decreased by a factor of 1.45 and if the motor driver is turned off, it is decreased by a factor of 1.26.

Conclusion on vibration analysis

In conclusion, the motor driver and the camera fans seem responsible for the motion of the PSF and this motion of the sources can be transmitted to the beam train thanks to different possible loose or unstable optics: the two DMs, the flat mirror replacing the SM, and the apodizer. All these optics are located before the FPM, which makes sense since when the coronagraph is on, the shape of the core of the PSF is changing. This study still gives a large panel of suspects and careful checking of the assembly and mounting screws for all testbed optomechanics did not result in any improvements in stability. Rather it appears that one or more mount designs themselves are insufficiently stable.

Mechanical improvement on the DM mount

As we showed in the previous part, DM1 seems to be the most likely source of the vibration effect, which can easily be explained since the current mount is really cantilevered. After this analysis, it was decided to re-design this mount to make it as stable as possible. The main constraints that have to be respected are the radius of curvature of the goniometer, the beam height (both of them have to correspond to the DM center) and the mount stability. We finally reach the global design presented in Fig. 4, where all the components are set above each other, following a vertical axis that makes the global structure really more stable than the current one.

SPECKLE NULLING SIMULATION TESTS

To study the instability phenomenon, we decided to focus on a well-known procedure: the Speckle Nulling method. Since Brown and Burrows 15 set the typical requirements in term of contrast for exoplanet detection, the theoretical feasibility of starlight subtraction in a so-called dark hole thanks to a DM has been proven. [START_REF] Malbet | High-Dynamic-Range Imaging Using a Deformable Mirror for Space Coronography[END_REF] This approach has been simplified with linearization of the equations to become the Speckle Nulling algorithm which has already several times been experimentally tested and has proven its efficiency and its robustness in broadband light. [START_REF] Trauger | Coronagraph contrast demonstrations with the high-contrast imaging testbed[END_REF][START_REF] Belikov | Toward 10 10 contrast for terrestrial exoplanet detection: demonstration of wavefront correction in a shaped-pupil coronagraph[END_REF] A first generalization of this method (Speckle Field Nulling), which is based on minimization of the speckle energy in all over the dark hole thanks to a Fast-Fourier-Transform-based (FFT) algorithm, has been developed. [START_REF] Bordé | High-Contrast Imaging from Space: Speckle Nulling in a Low-Aberration Regime[END_REF] In opposition to these estimation-free approaches, model-based techniques are developed and implemented. They separate estimation and control for a faster correction, such as the Self-Coherent Camera (SCC) [START_REF] Mazoyer | High-contrast imaging in polychromatic light with the self-coherent camera[END_REF] and the COronagraphic Focal-plane waveFront Estimation for Exoplanet detection (COFFEE) [START_REF] Paul | Simultaneous phase and amplitude retrieval with COFFEE: from theory to laboratory results[END_REF] for the estimation, or the Electric Field Conjugation (EFC) [START_REF] Riggs | Recursive starlight and bias estimation for high-contrast imaging with an extended Kalman filter[END_REF][START_REF] Give'on | Pair-wise, deformable mirror, image plane-based diversity electric field estimation for high contrast coronagraphy[END_REF] and the Stroke Minimization, [START_REF] Pueyo | Optimal dark hole generation via two deformable mirrors with stroke minimization[END_REF] including even multiple DMs to enable symmetric correction in the dark hole (both amplitude and phase aberration correction). In this section, we describe in the first part the Speckle Nulling code that is applied in simulation mode in the case of a circular aperture and the results we obtain. In the second part, we add random tip-tilt vibrations in the pupil plane of our simulated system and study the influence of this perturbation on the results.

Simulation without perturbation 4.1.1 Theory and code

We wrote a code that is split in different parts, to separate the simulation-only sections from the sections common to the simulation mode and to the on-testbed mode. The simulation-only sections are equivalent on the testbed to an image acquisition and therefore contain a simulated Lyot coronagraph similar to the HiCAT configuration. A visual description of the code is given in Fig. 5.

The Speckle Nulling method does not uniformly correct the image in the all dark hole. It typically focuses on the correction of the n brightest speckles of the dark hole. For example, if we only want to correct for one single speckle, the Speckle Nulling method selects the brightest speckle of the dark hole and computes its intensity I and its position (u x , u y ), according to the center of the PSF. Since one speckle in the field is equivalent to a sine function of the wavefront surface in the pupil plane, we are looking for an equation of the unknown error phase in the pupil plane such as:

Φ = A × sin(φ 0 + 2π × (f x x + f y y)) (1)
where A is the amplitude of the sine, φ 0 is the origin phase and f x and f y are the spatial frequencies of the phase and can easily be computed from the position (u x , u y ) of the speckle. Furthermore, to estimate the amplitude of the sine function, we add a calibration step before the correction. A known sine function is applied on the DM surface and we simulate the resulting image in the focal plane behind the coronagraph. This image has two symmetrical speckles, and we can get their intensity. According to this value and the amplitude of the sine command sent to the DM, we have a conversion factor from the DM command to the intensity of the image points.

A = C × A 0 × I I 0 (2)
where A is the calculated amplitude, A 0 is the amplitude of the calibration sine command, C is the conversion factor, I is the intensity of the speckle that has to be corrected, and I 0 is the intensity of the calibration speckle.

The main problem is finding the spatial phase φ 0 . As explained in Fig. 5, it is found by testing different possible spatial phases and selecting thanks to an interpolation the one that gives the best results on the corrected image.

Results

We apply the procedure described before on a simulated coronagraph with a Lyot Stop size equal to 0.99 times the pupil size. In this case, since we are working on simulated data, we fix the DM size to 136 actuators, so 4 times larger than our real DM. The error phase is set as shown on Fig. 6 on the top left image and we just have access to the top right image, which is the input of the routine.

We can also observe the results that we obtain after applying the Speckle Nulling code with the best correction to the images. As we can see, the correction was very efficient: in this case, the value of the speckle changed from 0.6898 to 0.01247 and on the corrected image, no speckle is visible. If we look at Fig. 7, the absence of speckle on the corrected image also clearly appears.

We also notice that on the simulated resulting phase in the pupil plane, issued from both the error phase and the correction phase, we can observe a non uniform phenomenon. This is due to the Moiré effect (frequency folding), which results from the difference between the spatial frequencies from the error phase and from the correction phase. Let's set the error phase has (f x , f y ) = (21, 21) as a spatial frequency and the frequency computed thanks to our Speckle Nulling code is (f x , f y ). The difference between the error phase frequencies and the correction frequencies is responsible from this Moiré effect and creates another periodic pattern in our pupil plane at a frequency (

|fx-f x | 2 , |fy-f y | 2
). This effect is negligible, since at the scale of the pupil plane, we just have a tiny fraction of this new period.

Simulation with perturbation

In this part, we add fake random vibration in the pupil plane before the coronagraph, that is equivalent of the DM plane.

Vibration simulation

Simulated vibratory perturbation can be applied on different components, such as the DM, the FPM or the Lyot Stop. Since on HiCAT the motion most likely comes from a component located before the FPM, we simulate vibrations in the DM plane. To do so, we generate a vector of n random tip values and a vector of n random tilt values, in a certain range. We can then compute the n corresponding wavefronts in the DM plane and we add them to our phase error. After that, we can simulate the n images behind the coronagraph that are used as an input of the procedure. We apply random vibration before summing the input images and also before applying the correction phases tested on the DM.

Except for these two extra steps, the code is applied similarly than in the previous case. We describe in the following part the results we obtain.

Results

As shown on Fig. 8, the vibration effect is clearly visible both before and after correction. In average, the correction is still doable, the vibration does not prevent the code from working. The main effects of the vibration are that we can obtain random results, sometimes really good and sometimes really bad, which is due to the fact that after correction, we take a single snapshot and the random tip-tilt can Right: one single snapshot after correction, with a similar random tip-tilt vibration. In both images, the tip and tilt aberrations are visible with the asymmetry of the diffraction rings that are produced by the Lyot coronagraph. As the correction is applied at a certain position, it does not necessarily correspond to the right speckle location since a tip-tilt is added. The correction is therefore less efficient than in the case we work without vibration. Figure 9. Azimuthal averaged intensity profiles of the simulated coronagraphic images in log scale for the noncorrected image (purple) and the corrected image (blue) that are displayed in Figure 8. The IWA and OWA of the controlled region are represented with vertical lines. A shift of the diffraction rings is observed between the two cases. The attenuation of the speckle at its location is less important than the reduction observed in the absence of vibration in Figure 7.

be small or important. Furthermore, the translation of the pattern appears on the azimuthal plot. Indeed, we observe both a translation of the rings and a smoothing of the values, as shown on Fig. 9.

The image 10 indicates the performance depending on the vibration amplitude maximum range. For every maximum range, we apply the code 5 times and obtain a set of 5 speckle values after correction. We then look at the minimum obtained value, which corresponds to the best performance, the average value and the maximum value (worst correction). As shown on the plot 10, the higher the instability amplitude is, the more random the results are, even if in every case we can still obtain very good results (the minimum value is quite constant).

Figure 10. Intensity of the speckle after correction as a function of the maximum vibration amplitude. For every maximum range, we apply the code 5 times and obtain a set of 5 speckle values after correction. From this analysis, we retrieve the minimum obtained value, which corresponds to the best performance (yellow curve), the average value (blue curve) and the maximum value that corresponds to the worst correction (pink curve).

CONCLUSIONS

The HiCAT testbed will enable high-contrast imaging studies for telescopes with segmented apertures, spiders and/or central obstruction, such as WFIRST or HDST. It is designed as a coronagraph completed with two deformable mirrors to perform wavefront control and address both phase and amplitude aberrations. This year, we implemented first wavefront control codes that led us to detect an instability issue with drifts that were observed on the PSF. This issue needs to be solved to pursue our wavefront control studies. Thanks to the PSF motion analysis, we obtained strong suspects that could cause this instability. Further studies will be done thanks to an accelerometer to confirm our initial guess. In the meanwhile, a new DM mount has been designed and should provide an important improvement in terms of testbed stability. In parallel, we coded a preliminary version of a Speckle Nulling code that gives good results in simulation mode, even in presence of pupil plane instability. Without vibration, a given speckle is fully removed after correction. In presence of vibration, we obtain a reduction of the speckle intensity in almost every case. This code will further be applied on the testbed to correct for the speckles in the chosen dark zone. Our code will also be adapted to the case of a segmented mirror to study the effects of segments that vibrate independently on contrast performance.

Figure 1 .

 1 Figure1. HiCAT testbed design done with the software Solidworks, the beam is exported from Zemax. The telescope is simulated by a pupil mask, the segmented mirror and off-axis parabolas. The segmented mirror is conjugated with the pupil mask to form a segmented pupil with central obstruction and spider struts. The off-axis parabolas set the telescope aperture. The wavefront control is done with two deformable mirrors. The coronagraph is composed of an apodizer, a focal plane mask and a Lyot stop.

Figure 2 .

 2 Figure2. Analysis of the PSF motion when HiCAT is at rest (all the hardware is on, the enclosure is closed and the camera is at the end of the testbed). Left: distance in pixels of the PSF from the mean position as a function of time. Center: direction of the motion on the detector as a function of time. No privileged direction appears here. Right: Spectrum of the motion. The motion is obviously a low frequency movement, there is no high frequency peak.

Figure 3 .

 3 Figure 3. Temporal amplitude (top) and spectra (bottom) of the PSF motion obtained from a video acquisition when tapping three optical components. Left: the flat mirror replacing DM2, middle: the flat mirror replacing the apodizer, right: the flat mirror replacing the Iris-AO segmented mirror. The red lines correspond to the maximum value obtained without tapping.

Figure 4 .

 4 Figure 4. Solidworks design of the mount. Both the center of curvature of the goniometer and the beam height have an origin that corresponds to the DM center. The components are set one above each other, which guarantees better stability than the current cantilevered mount design.

Figure 5 .

 5 Figure 5. Diagram of our Speckle Nulling code. The code has been written to work in two different modes: with simulated data or experimental images on the testbed. Simulations 1 and 2 are dedicated to the simulation mode only while Codes 1, 2 and 3 are common to the simulation mode and the on-testbed mode. Simulation 1 sets up all the useful input for the rest of the procedure: a fake flat command of the DM, an "unknown" phase error (here: a sine function), background frames and image plane frames obtained with a Lyot coronagraph. Code 1 computes the characteristics of the selected speckle (spatial frequency and amplitude) and takes the average background-subtracted image. Code 2 generates n commands to the DM to probe different spatial phases of the selected speckle. Simulation 2 computes the coronagraphic images in the presence of unknown phase for each of the n commands applied on the DM. Finally, Code 3 analyses the coronagraphic images to determine the phase and the corresponding DM command that allow us to obtain the best correction of the selected speckle.

Figure 6 .

 6 Figure 6. First results obtained with the Speckle Nulling code with simulated images before and after correction of a pair of speckles (top and bottom). Left: Images in the relayed pupil plane of the coronagraph before application of the Lyot stop. Right: Coronagraphic images in the final image plane.

Figure 7 .

 7 Figure 7. Azimuthal averaged intensity profiles of the coronagraphic images in log scale for the initial image (purple) and the corrected image (blue) that are displayed in Figure 6. The IWA and OWA of the controlled region are represented with vertical lines. The attenuation of the intensity is observed at the speckle location ( 44 pixels).

Figure 8 .

 8 Figure 8. Coronagraphic images without and with correction, in presence of vibration. Left: Averaged image of 10 frames with random tip and tilt within the amplitude range [-0.01λ, +0.01λ].Right: one single snapshot after correction, with a similar random tip-tilt vibration. In both images, the tip and tilt aberrations are visible with the asymmetry of the diffraction rings that are produced by the Lyot coronagraph. As the correction is applied at a certain position, it does not necessarily correspond to the right speckle location since a tip-tilt is added. The correction is therefore less efficient than in the case we work without vibration.