domain: spi phys

text:
Daniel Destarac
email: daniel.destarac@onera.fr

Gérald Carrier
email: gerald.carrier@onera.fr

George R Anderson
email: george.anderson@nasa.gov

Siva Nadarajah
email: siva.nadarajah@mcgill.ca

Daniel J Poole
email: d.j.poole@bristol.ac.uk

John C Vassberg
email: john.c.vassberg@boeing.com

David W Zingg

An Example of a Pitfall in Aerodynamic Shape Optimization

AIAA) launched an optimization discussion group, ADO-DG, in 2013. One of the four benchmark test cases is based on the NACA0012 airfoil and solutions of the Euler equations with prescribed objective function and geometric constaints. Volunteer participants were invited by ONERA to submit their optimized shapes for uniform pressure drag assessment through grid convergence. The assessment required three phases. Initial investigations with a simple protocol proved insufficient and raised some questions. This led to further investigations which showed that for most optimized airfoils there exists a lower drag branch and a higher drag branch in the drag (Mach number) diagram. Hysteresis was observed when performing a downward Mach number sweep and an upward Mach number sweep. For several airfoils, the design Mach number falls within the overlap range of the two branches, which indicates that two solutions may exist for the same airfoil. Two kinds of supersonic flow structure can be distinguished. One kind produces a lower drag branch, the other a higher drag branch. The jump from one branch to the other involves non symmetrical converged solutions.

These further investigations required a final assessment, where a protocol was instituted to ensure that the lower drag branch was captured by taking advantage of the hysteresis. However, nonuniqueness may confound an optimization algorithm, as the same set of design variables can give two different objective function values.

I. Introduction

AIAA launched in 2013 the Aero-Design Optimization Discussion Group ADO-DG with prescribed common test cases. One of the test cases is based on the NACA0012 airfoil with solutions of the Euler equations. The objective function to be minimized is the pressure drag coefficient CD p and the constraint is that at any chordwise location, thickness must be greater than, or equal to, that of the initial airfoil. The aerodynamic conditions are for a single design point at M ∞ = 0.85, α = 0 • . Volunteer participants were invited by ONERA to submit their optimized shapes as a list of point coordinates for uniform pressure drag assessment through grid convergence. Seven shapes from six participants, ONERA, NASA/Stanford University, McGill University, University of Bristol, The Boeing Company and University of Toronto, were submitted.

They are documented in references [START_REF] Carrier | Gradient-Based Aerodynamic Optimization with the elsA Software[END_REF][START_REF] Bisson | Adjoint-Based Aerodynamic Optimization of Benchmark Problems[END_REF][START_REF] Bisson | Adjoint-Based Aerodynamic Optimization of Benchmark Problems[END_REF][START_REF] Hicken | Aerodynamic Optimization Algorithm with Integrated Geometry Parameterization and Mesh Movement[END_REF][START_REF] Telidetzki | Application of Jetstream to a Suite of Aerodynamic Shape Optimization Problems[END_REF][START_REF] Lee | Aerodynamic Shape Optimization of Benchmark Problems Using Jetstream[END_REF][START_REF] Poole | Control Point-Based Aerodynamic Shape Optimization Applied to AIAA ADODG Test Cases," to Two-Dimensional Drag Minimization[END_REF][START_REF] Ledoux | A Study Based on the AIAA Aerodynamic Design Optimization Workshop Test Cases[END_REF][START_REF] Anderson | Aerodynamic Shape Optimization Benchmarks with Error Control and Automatic Parameterization[END_REF][START_REF] Méheut | Gradient-Based Single and Multi-Point Aerodynamic Optimizations with the elsA Software[END_REF]. They were obtained with very different parameterization types, numbers of design variables, optimization algorithms and flow solvers. It seemed in line with the Discussion Group spirit that one partner carry out an assessment process, identical for all the optimized airfoils.

In Section 2 baseline and optimized airfoils are presented. In Section 3 the planned assessment protocol is described. In Section 4 it is shown that for such unusual airfoils as those investigated here, the planned protocol is too simple. In Section 5, further investigations provide clearer insight into the phenomenas. This insight leads to the final assessment presented in Section 6. Concluding remarks are given in Section 7.

II. Baseline and Optimized Shapes

The test case [11] is based on work by Vassberg et al. [START_REF] Vassberg | A Systematic Study on the Impact of Dimensionality for a Two-Dimensional Aerodynamic Optimization Model Problem[END_REF]. However, here trailing edge closure at x = 1 is ensured through a modification of the polynomial definition of the NACA0012 x 4 coefficient. The modified airfoil equation is y = ±0.6 (0.2969 √ x -0.1260 x -0.3516 x 2 + 0.2843 x 3 -0.1036 x 4 )

Information about the optimized airfoils and the approaches used is contained in table 1.

The upper lines of the baseline and optimized airfoils are plotted in figures 1 and 2. All optimized airfoils sat- isfy the thickness constraint, thickness greater or equal to the baseline airfoil thickness at all chordwise locations. All optimized airfoils have common features: a blunt leading edge, a "flat top", a highly curved rear region and a large trailing edge included angle, which were also features of the optimized airfoils in [START_REF] Vassberg | A Systematic Study on the Impact of Dimensionality for a Two-Dimensional Aerodynamic Optimization Model Problem[END_REF]. However, as expected given the variety of the approaches (table 1), significant differences appear between their shapes.

III. Planned Assessment Protocol

The grids are based on some of the NACA0012 airfoil grids of Vassberg and Jameson [START_REF] Vassberg | In Pursuit of Grid Convergence for Two-Dimensional Euler Solutions[END_REF]. These highly regular grids with quasi-orthogonal cells have an O-type topology with identical numbers of cells in the two directions. If nc is this number, levels nc = 256, 512, 1024 are used in this study. The grid for level nc = 256 for airfoil 1 is illustrated in figure 4. The metric properties of these grids were designed to ensure a quasi linear behavior of pressure drag with 1/nc 2 for usual second order finite volume formulations. That was the case in [START_REF] Vassberg | In Pursuit of Grid Convergence for Two-Dimensional Euler Solutions[END_REF], where extrapolation to zero mesh size was thus straightforward. It remained so for the various optimized airfoils produced by ONERA and compared in [START_REF] Carrier | Gradient-Based Aerodynamic Optimization with the elsA Software[END_REF].

The optimized airfoils, requested to number at least 129 points, are transformed into series of 2048 points through spline interpolation. The NACA0012 nc = 2048 grid is then projected onto their shape. The resulting grid is coarsened to levels nc = 1024, 512 and 256 for the resolution of the Euler equations.

The Euler solutions are computed with the ONERA elsA software [START_REF] Cambier | Status of the elsA Software for Flow Simulation and Multi-Disciplinary Applications[END_REF][START_REF] Cambier | The Onera elsA CFD Software: Input from Research and Feedback from Industry[END_REF]. As in the assessment computations of [START_REF] Carrier | Gradient-Based Aerodynamic Optimization with the elsA Software[END_REF], the centred scheme of Jameson [START_REF] Jameson | Numerical Solution of the Euler Equations by Finite Volume Methods Using Runge-Kutta Time-Stepping Schemes[END_REF] with the JST artificial dissipation is used, dissipation coefficients being set to κ2 = 0.5, κ4 = 0.008. To facilitate convergence, a wall slip condition instead of a matching condition is applied to the grid lines joining the trailing edge and the downstream grid boundary, see figure 4. Based on previous experience, the default number of iterations for the present study was set to 5000, with 3 multigrid levels.

The planned protocol was, for each airfoil, to perform three computations, nc = 256, 512, 1024, at the unique design conditions, M ∞ = 0.85, α = 0 • , and to extrapolate to zero mesh size through the procedure described in [START_REF] Vassberg | In Pursuit of Grid Convergence for Two-Dimensional Euler Solutions[END_REF].

Things did not turn out as straightforwardly as anticipated and led to further investigations.

IV. Initial Investigations

The planned protocol was first systematically applied to the seven airfoils. Anomalies appeared for some of them. At this stage, it is enough to consider airfoils 1 and 2.

Convergence of the residual L ∞ norm of variable ρ is plotted in figures 5 and 6. For airfoil 1, convergence is as expected from common experience, reaching machine zero level and requiring more iterations with a finer grid level than with a coarser (figure 5). Such is not the case for airfoil 2 (figure 6), where machine zero level is reached with the coarsest and the finest grid but not with the medium grid. Pressure drag computation was still performed, though with dubious results. Indeed a striking anomaly affects the behavior of airfoil 2 pressure drag through grid refinement, plotted in figure 7. Pressure drag of airfoil 1 decreases with grid refinement because spurious drag [START_REF] Destarac | Far-Field / Near-Field Drag Balance and Applications of Drag Extraction in CFD[END_REF] (mostly at the blunt leading edge) decreases, while the flow physics (shock waves) are not disturbed. The opposite behavior with airfoil 2, added to the convergence anomaly, leads us to suspect some flow topology change with grid refinement.

Chordwise pressure coefficient distributions on airfoil 2 for the three grid refinement levels are plotted in figure 8.

The topology change is from a unique supersonic/subsonic shock trace with the coarsest grid (thick line in the figure )   to an initial supersonic/subsonic shock trace followed by a strong expansion leading to a second supersonic/subsonic shock trace with the medium and fine grid.

This calls for several remarks: 1) the double shock pattern was not found by the designers of airfoil 2 who checked it with several solvers [START_REF] Bisson | Adjoint-Based Aerodynamic Optimization of Benchmark Problems[END_REF] and found consistently much lower drag coefficient levels (around 40 drag counts, not 90 as in figure 7); 2) the difficult convergence with the medium grid must come from the fact that the flow is still close to the topology change; 3) the double shock will produce much more wave drag than the unique shock (in the present case), which accounts for the anomaly of figure 7, which shows much higher drag with fine grid than with coarse grid.

As pointed out by Vassberg et al. in [START_REF] Vassberg | A Systematic Study on the Impact of Dimensionality for a Two-Dimensional Aerodynamic Optimization Model Problem[END_REF], to understand drag production with airfoils of the kind dealt with here, "one must inspect the flowfield, not just the properties on the airfoil surface". This will be done extensively in a subsequent section entitled Final assessment. Here such inspection is limited to the finest grid solutions for airfoils 1 and 2, figures 9 and 10 displaying supersonic iso-Mach number line patterns. Airfoil 1 shows a single supersonic area (for a half-airfoil), airfoil 2 a double supersonic area. Airfoil 1 produces a supersonic/subsonic shock at a distance in the field and, near the airfoil, a more complex structure consisting of a supersonic/supersonic shock followed by a short supersonic expansion and ending with a supersonic/subsonic shock. Airfoil 2 produces an extended supersonic/subsonic shock in the field and a subsonic-to-supersonic expansion leading to a second less extended supersonic/subsonic shock.

The first type of solution will be called solution of the first kind in this paper and the second type solution of the second kind.

V. Further Investigations

Airfoil 2 supports a solution of the first kind in [START_REF] Bisson | Adjoint-Based Aerodynamic Optimization of Benchmark Problems[END_REF] and with the coarsest grid of the present study, but a solution of the second kind with the finer grids. At this stage, an hypothesis could not be neglected, that of "non-unique solutions of a discrete model of the Euler equations" as first evidenced by Jameson [START_REF] Jameson | Airfoils Admitting Non-Unique Solutions to the Euler Equations[END_REF], then by Hafez and Guo [START_REF] Hafez | Nonuniqueness of Transonic Flows[END_REF]. In [START_REF] Jameson | Further Studies of Airfoils Supporting Non-Unique Solutions in Transonic Flow[END_REF],

Jameson et al. produced other non-unique solutions for airfoils of the same type as those investigated in the present paper. Ou et al. [START_REF] Ou | Airfoils Supporting Non-unique Transonic Solutions for Unsteady Viscous Flows[END_REF] confirmed multiple solutions for airfoils from [START_REF] Jameson | Further Studies of Airfoils Supporting Non-Unique Solutions in Transonic Flow[END_REF] in unsteady viscous flow. The same authors showed that the problem is not limited to 2D airfoils but may arise for wings as well [START_REF] Ou | Studies of Wings Supporting Non-unique Solutions in Transonic Flows[END_REF].

Similar shock structures to those encountered in the present study are also discussed by Kuz'min and Ivanova in [START_REF] Kuz'min | The structural Instability of Transonic Flow Associated with Amalgamation/Splitting of Supersonic Regions[END_REF]: "singular free-stream Mach numbers, which trigger the splitting/amalgamation of local supersonic regions".

A physical interpretation of this phenomenon given by these authors is that in non-isentropic steady flow, the downstream boundary of a supersonic area cannot physically be in contact with the upstream boundary of another supersonic area. So, when, through some infinitesimal perturbation, a unique supersonic area is split into two areas, the process must be discontinuous.

In these further investigations, it was chosen to take Mach number as the perturbation factor. Downward Mach number sweeps were computed with the coarsest nc = 256 grid level for all airfoils, around the design Mach number first kind to a solution of the second kind through the splitting of a single supersonic area into two as described in [START_REF] Kuz'min | The structural Instability of Transonic Flow Associated with Amalgamation/Splitting of Supersonic Regions[END_REF],

M ∞ = 0.
with a correlative discontinuous drag increase. This occurs just below the design point Mach number M ∞ = 0.85. It can be seen that airfoil 5 behaves differently from the other airfoils. The reason is that, for this airfoil, transition from a solution of the first kind to a solution of the second kind does not happen through the splitting of a single supersonic area but through the appearance of a second, separate, continuously growing, supersonic area, with no contact with the primary one in the investigated range.

From these coarse grid results, three airfoils were selected for fine grid checking (nc = 1024), and performing It should be noted that a similar behavior was found for one of the present airfoils with an entirely different flow solver [START_REF] Lee | Aerodynamic Shape Optimization of Benchmark Problems Using Jetstream[END_REF]: it is not an artefact of the solver used in this study.

The shock pattern of solutions of the first kind for airfoils 2 and 7 shows a long supersonic-subsonic wave away from the airfoil surface, broken down into a short supersonic-supersonic wave and a short supersonic-subsonic wave near the airfoil surface. This might be the feature which accounts at once for very low drag and instability. In solutions of the second kind, the pattern is simpler: two supersonic-subsonic waves separated by an expansion.

VI. Final Assessment

For airfoils 2 and 7 the existence of a lower drag branch over the design Mach number and of a higher drag branch below, with some overlap providing double solutions, has thus been evidenced. The lower drag branch consists of solutions of the first kind while the upper drag branch consists of solutions of the second kind. Figures 11 and12 suggest that it might also be the case for airfoils 1, 3, 4, 6. So a modified assessment protocol taking advantage of the downward sweep hysteresis was devised: 2000 iterations at M ∞ = 0.8505 + 2000 iterations at M ∞ = 0.85025 + 5000 iterations at M ∞ = 0.85. This procedure led computations of all airfoils susceptible of a double branch (1, 2, 3, 4, 6, 7) to lower drag branch solutions of the first kind, see figures 18, 19, 20, 21, 23, 24. In the case of airfoil 5, which does not have a branch at the design point, the solution shows a single supersonic area (for a half-airfoil), figure 22, but the iso-Mach number line pattern is not strictly of the first kind. A close look at figure 22 does show a short supersonic/supersonic wave, but this wave is detached from the airfoil.

Computed pressure drag coefficient values are gathered in table 2 and figure 25. With the modified assessment protocol, for all airfoils drag decreases as the grid is being refined, compare figure 25 to figure 7. Extrapolation to zero mesh size was attempted following the method of [START_REF] Vassberg | In Pursuit of Grid Convergence for Two-Dimensional Euler Solutions[END_REF]. The measured order of accuracy with this method varies between 1.5 and 2.2 depending on the airfoil considered. But considering the coarseness of level nc = 256 and the imperfect linear behavior for some airfoils, extrapolated data were not retained.

Different grid dependency slopes in figure 25 reflect different levels of spurious drag. More or less blunt leading edges will induce more or less spurious drag because spurious drag at the leading edge is influenced by pressure gradients. This can be illustrated by comparing the curvature of airfoils 3 and 6 in figure 26. Curvature is computed with the reduced curvilinear abscissa parameter, 0 at the lower side trailing edge, 0.5 at the leading edge, 1 at the upper side trailing edge. Airfoil 6 is much blunter than airfoil 3. Consistently, the grid dependency slope in figure 25 is much higher for airfoil 6 than for airfoil 3. A normally rounded leading edge such as the NACA0012's produces weaker pressure gradients than the blunt leading edges of the airfoils investigated here. (With an inviscid flow model, a squared leading edge would induce infinite pressure gradients.) The much larger amounts of spurious drag with the optimized shapes than with the NACA0012 is obvious in table 2.

If the differences in grid dependency slope arise from numerical effects (spurious drag), those in finest grid drag level arise mostly from physical phenomena (wave drag). All airfoils but one support a shock structure of the first kind, airfoils 1, 2, 3, 4, 6, 7. However there are differences between these structures (figures [START_REF] Destarac | Far-Field / Near-Field Drag Balance and Applications of Drag Extraction in CFD[END_REF][START_REF] Jameson | Airfoils Admitting Non-Unique Solutions to the Euler Equations[END_REF][START_REF] Hafez | Nonuniqueness of Transonic Flows[END_REF][START_REF] Jameson | Further Studies of Airfoils Supporting Non-Unique Solutions in Transonic Flow[END_REF][START_REF] Ou | Studies of Wings Supporting Non-unique Solutions in Transonic Flows[END_REF][START_REF] Kuz'min | The structural Instability of Transonic Flow Associated with Amalgamation/Splitting of Supersonic Regions[END_REF].

An attempt at identifying a relation between drag and the length of the supersonic / supersonic wave is presented in figure 27 with an exponential fit at the data points. It is acknowledged that the length measurement involves a part of subjectivity, but the ambition is only to identify this supersonic wave length as a major factor in wave drag production in solutions of the first kind. It is obvious that there must be other shock features that contribute to the wave drag variation, such as shock angle and curvature.

VII. Conclusion

One of the test-cases of the AIAA Aero-Design Optimization Discussion Group is based on the NACA0012 airfoil with solutions of the Euler equations. The purpose of the present study was to provide a pressure drag assessmant of the optimized shapes by several participants through grid convergence. It had been anticipated that this study would be straightforward and dull. Fortunately, initial investigations produced anomalies and raised questions. Further investigations, much inspired by Jameson's 1991 paper on non-unique numerical solutions of the Euler equations and Jameson et al. sequel of 2012, showed that for most of the airfoils investigated there exists a lower drag branch and a higher drag branch in the drag coefficient vs. Mach number diagram. Hysteresis when performing a downward Mach number sweep and an upward Mach number sweep was such that, for several airfoils, the design point was inside the overlap range of the two branches, which means a double solution at this point. A low drag converged solution or a high drag converged solution could be found depending on the history of the computation. Dependency on computation history is rather unwelcome in CFD but cannot be ignored. Two typical structures of supersonic flow were distinguished. One consists in a single supersonic flow area containing a supersonic/supersonic wave starting from the airfoil surface. The other involves two separate supersonic flow areas without a supersonic/supersonic wave. The jump from one branch to the other is due to flow topology change from one kind to the other. The first kind produces a lower drag branch, the second kind a higher drag branch. The topology change from one kind to the other may involve non symmetrical converged solutions, of the first kind on one side of the airfoil and of the second kind on the other.

For the final assessment, a procedure that takes advantage of the hysteresis was used to ensure capture of the lower drag branch at the design point. Reliable drag dependency on grid refinement was thus obtained. A relation between the drag efficiency and the length of the supersonic/supersonic wave has been proposed. All airfoils but one admit two branches, the optimum being a solution of the first kind. The airfoil behaving differently is of neither kind. It supports a single supersonic area but the supersonic wave is separated from the airfoil. Its drag level is among the highest, but its behavior in a Mach sweep is continuous, as if excellence in single-point drag minimization in highly transonic flow involved the drawback of instability in the vicinity of the optimum. The conclusion of Jameson's 1991 paper contains the warning that optimization does not necessarily lead to a good design. As regards nonunique solutions, which have been identified by other authors, although there is no reason why they might not occur with shapes not derived from optimization, the present paper indicates a strong correlation between optimality and nonuniqueness: six of the seven optimized airfoils support nonunique solutions.

For the airfoils apt to support multiple solutions, the optimization algorithms used by the various participants always led to the lower drag branch, which is satisfactory for the ADO-DG purpose. However, nonuniqueness may have negative effect on convergence of an optimization algorithm (as the same set of design variables can give two different objective function values). The present paper, concentrating on numerical flow physics investigations, does not contain a critical discussion of the different optimization methods or parameterizations. From the data produced here, a follow-up paper dedicated to such a discussion might be of interest. However the present paper already provides an overview of the status of the optimization results obtained by the different teams having worked on the ADO-DG test case # 1 during the two first years since its creation. It will also, hopefully, encourage new participants to investigate this optimization test case and try to achieve designs even closer to a shock-free airfoil.

Figures

  an upward Mach number sweep after the downward sweep. The objectives were to search for hysteresis behavior which would account for multiple solutions, and to investigate probable transient non-symmetrical solutions (one appears in figure12for airfoil 7). The chosen airfoils are airfoil 2 because of suspicion of multiple solutions, airfoil 5 for its unique continuous behavior and airfoil 7 the lowest drag airfoil at design conditions. The width of the explored Mach number range is narrower than with coarse grids, between M ∞ = .8490 and M ∞ = .8505. The basic step, ∆M ∞ = ±0.0001, has been empirically refined around discontinuities for airfoil 7 in search of non-symmetrical solutions.CD p (M ∞ ) curves are plotted in figure13. For airfoils 2 and 7, this figure shows a lower drag branch and a higher drag branch. During the downward Mach number sweep, the solution jumps from the lower branch to the higher branch, while in the upward sweep it jumps back to the lower drag branch. Hysteresis makes the two jump Mach numbers different. So, in a narrow Mach number band, double solutions may be obtained. The design Mach number, M ∞ = 0.85 is inside this band. The lower drag branch consists of solutions of the first kind, the upper drag branch of solutions of the second kind. With airfoil 5 a unique branch is found, with continuous transition from the first to the second kind through the development of a downstream supersonic area separated from the primary one. For airfoil 7, figures 14, 15, 16 illustrate the two-phase topology change from the first kind to the second kind through a non-symmetrical solution, first kind on the upper side, second kind on the lower side. Convergence of these three computations is plotted in figure17. It is slow immediately before the topology change, fast in the nonsymmetrical case and immediately after the topology change. In the upward sweep, the non-symmetrical solution (second kind on the upper side, first kind on the lower side) is less transient than in the downward sweep.

Fig. 1 Fig. 2 Fig. 3 Fig. 4 Fig. 5 Fig. 6 Fig. 7 Fig. 8 Fig. 9 Fig. 10 Fig. 11 Fig. 12 7 MFig. 13 FineFig. 14 Fig. 15 Fig. 16 Fig. 17 Fig. 18 Fig. 19 Fig. 20 Fig. 21 Fig. 22 Fig. 23 Fig. 24 Fig. 25 Fig. 27

 12345678910111271314151617181920212223242527 Fig.1Baseline and optimized shapes[START_REF] Carrier | Gradient-Based Aerodynamic Optimization with the elsA Software[END_REF]

Table 1

 1 Information about optimized airfoils

	identification	design	parameter	design variables	optimization	optimization
			-ization type		algorithm	grid size
						(cell numbers)
	airfoil 1 airfoil 2	ONERA McGill	Bézier B-splines	96 → 127 16	SLSQP → DOT SNOPT [13]	65,536 98,000
		University
	airfoil 3	University	B-splines	18	SNOPT [13]	192,000
		of Toronto
	airfoil 4 airfoil 5	University of Bristol Boeing	Symmetric SVD modes Free Surface	8 → 15 799	GSA (8 modes) → FSQP (15 modes) Smoothed	131,072 1,310,720
			(Symmetric)		Steepest Descent
	airfoil 6	NASA / Stanford	radial basis functions	7 → 15 → 31	SNOPT [13]	122,000
		University
	airfoil 7	ONERA	Bézier	96	airfoil 1
					→ SLSQP	65,536

  from the converged state at M ∞ + ∆M ∞ . The prescribed number of iterations was 3000. Most computations reached machine zero. When such was not the case, solutions with uncertainty |∆CD p | ≥ .01 drag count were rejected, hence the gaps seen in figures 11 and 12. The purpose of these tests was looking for solutions of the first and second kind and the topology change from one to the other. If an airfoil may support both kinds of solution, high Mach numbers are likely to produce the first kind, with a shock located close to the trailing edge. Indeed, for all airfoils, a solution of the first kind is obtained with the maximum Mach number M ∞ = 0.8505. When the Mach number is lowered, a solution of the second type is always obtained, but through two different mechanisms, one discontinuous, one continuous, see figures 11 and 12. All airfoils but airfoil 5 undergo a discontinuous jump from a solution of the

85 with a step ∆M ∞ = -0.0001. Except for the first computation, at M ∞ = 0.85, all subsequent calculations were started

Table 2

 2 Grid convergence of pressure drag coefficient (expressed in drag counts)

	mesh size	CD p	CD p	CD p	CD p	CD p	CD p	CD p	CD p
		NACA0012 airfoil 1 airfoil 2 airfoil 3 airfoil 4 airfoil 5 airfoil 6 airfoil 7
	nc = 256	473.3	58.3	69.0	54.0	106.3	99.6	84.9	52.3
	nc = 512	471.6	40.8	40.6	44.5	89.8	87.7	55.1	32.5
	nc = 1024	471.2	36.7	31.8	42.4	86.1	85.0	44.7	27.9


domain: spi phys math info

text:
Ioannis Sarras

Hélène Piet-Lahanier

Julien Marzat

Sylvain Bertrand

Emmanuel Nuño

Comments on 'Distributed observers design for leader-following control of multi-agent networks

Keywords: Multi-agent system, Active leader, Distributed control, Distributed observer, Common Lyapunov function

In this note it is pointed out that the proof of Theorem 1 of [1] is incorrect for the theorem as stated. A proof for a modified version of the theorem is also given 1 .

The seminal contribution of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF] was the extension of the conventional observer design to a distributed observers' design and in particular for a scenario of a multi-agent system, where an active leader with unknown velocity has to be followed. This work has been extremely popular and, as a result, has been cited by an increasing number of publications on multi-agent systems. However, the proof of Theorem 1 of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF] is incorrect for the theorem as stated. In this note we give a proof for a modified version of the theorem. This yields a different set of sufficient conditions which establish the stability proof and highlight the restrictive conditions on the acceptable topologies. In what follows we adopt the same notation as in the original work [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF].

The dynamic models for the leader and the agents consist of double integrators as in equations ( 1)-( 2) of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF]. Then, This work was supported by the project PR-GUIMAUVE of ONERA. The last author is supported by the Mexican CONACyT basic scientific research grant CB-282807. This paper was not presented at any IFAC meeting. Corresponding author I. Sarras. Tel. +33-1-80-386665.

Email addresses: ioannis.sarras@onera.fr (Ioannis Sarras), helene.piet-lahanier@onera.fr (Hélène Piet-Lahanier), julien.marzat@onera.fr (Julien Marzat), sylvain.bertrand@onera.fr (Sylvain Bertrand), emmanuel.nuno@cucei.udg.mx (Emmanuel Nuño). 1 While the paper was under review, we became aware that another way to address the error in [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF] is to simply replace 1/k with 1/k 2 in equation (4) of that paper, as was done in [START_REF] Wang | Distributed observers for tracking a moving target by cooperative multiple agents with time delays[END_REF]. With such change the proof of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF] is valid. We now also understand that the simulations in [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF] were carried with the corrected 1/k 2 in equation (4) of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF].

the authors proposed the neighbor-based control for agent i, given in (3) of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF], to track the leader while also obtaining an estimate of the (unavailable) leader velocity by the observer in equation ( 4) of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF]. We first note that the expression of the closed-loop dynamics in equation ( 6) in the proof of Theorem 1 of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF] is incorrect since the k coefficient in the last term of η was missing. This puts into question the correctness of the statement of Theorem 1 in [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF], and to this end we propose an alternative stabilization theorem the for coordination rule (3)-(4) of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF]. Theorem 1 Consider the closed-loop multi-agent system. In each time interval [t i , t i+1 ), if the entire graph is connected and satisfies

λ > 8 9 λ, (1)
then, there exist gains k > 0 and l > 0 satisfying the conditions

k > max 2, 1 2 1 + 1 + 8λ 2 9λ 2 -8λ 2 (2)
l 1 ≤ l ≤ l 2 , (3)
with PROOF. We will provide the stability proof for one time interval as in the original work. Following the notation in [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF] we express the closed-loop dynamics in the compact form of equation ( 6) of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF] but with the matrix

l 1,2 := 12k 2 (k -1)λ ± √ ∆ 2(2k -1) 2 λ 2 ∆ := -16k 3 (k -1) 8λ 2 k 2 -(9λ 2 + 8λ 2 )k + 2λ
F p :=     0 I 0 -lH p -kI kI -l k H p 0 0     . (4)
Now, define the candidate Lyapunov function V as in equation (8), of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF], that is positive-definite for k > 1. The timederivative of V along z-trajectories gives (9) of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF] but with

Q p :=     lH p l 2k (2k -1)H p -kI l 2k (2k -1)H p 2(k -1)I -kI -kI -kI kI     .
To check positive definiteness of Q p , and ensure global exponential stability of the origin, we use the Schur's lemma with the matrix definitions

A := lH p l 2k (2k -1)H p l 2k (2k -1)H p 2(k -1)I , E := -kI -kI
, C := kI to obtain the inequalities

D := A -EC -1 E T = lH p -kI l(2k-1) 2k H p -kI l(2k-1) 2k H p -kI (k -2)I =: D 1 D 2 D 2 D 3 > 0, k > 0.
While the second one is trivially satisfied, for the positive-definiteness of matrix D we invoke, again, the Schur's lemma to obtain the conditions

D 3 > 0 ⇐⇒ k > 2 and D 1 -D 2 D -1 3 D T 2 > 0 ⇐⇒ (2k -1) 2 l 2 H 2 p -12k 2 (k -1)lH p + 8k 3 (k -1)I < 0.

In addition, using the bounds on the matrix

H p yields (2k -1) 2 λ 2 l 2 -12k 2 (k -1)λl + 8k 3 (k -1) < 0.
This last condition can be seen as a polynomial of second degree in l, expressed as f (l) := αl 2 + βl + γ, α > 0, β < 0, γ > 0, that needs to be negative, i.e. f (l) < 0. We will now derive the general conditions for which set of values of k, l this is true. We calculate first the discriminant

∆ = 16k 3 (k -1)(9λ 2 -8λ 2 ) k 2 -k -2λ 2 9λ 2 -8λ 2
, which needs to be positive to provide us with acceptable real solutions. For ∆ > 0, and taking into account also that k > 2, we have two cases.

(1) 9λ 2 -8λ 2 < 0: This condition is always satisfied since

λ 2 < 8 9 λ 2 < λ 2 .
For positivity of ∆ the following condition should hold g(k

) := k 2 -k -2λ 2 9λ 2 -8λ 2 < 0. Calculating the discriminant for g(k) = 0 yields ∆ k = 1 9λ 2 -8λ 2 9λ 2 that
can never be positive since 9λ 2 -8λ 2 < 0.

(2) 9λ 2 -8λ 2 > 0: In this case, the topologies are constrained to the configurations where λ 2 > 8 9 λ 2 . Now, for positivity of ∆ one requires g(k

) := k 2 -k -2λ 2 9λ 2 -8λ 2 > 0. The discriminant for g(k) = 0 reads ∆ k = 1 9λ 2 -8λ
2 9λ 2 and is always positive. We can deduce that g(k) = 0 has two roots, one negative k 1 and one positive k

2 , with |k 1 | < k 2 , with k 2 explicitly defined as k 2 := 1 2 1 + 1 + 8λ 2 9λ 2 -8λ 2 .
Hence, we obtain that, for g(k) > 0, k should satisfy

the condition k > 1 2 1 + 1 + 8λ 2 9λ 2 -8λ 2
. Now for the polynomial f (l), due to the fact that its coefficients are such that α > 0, β < 0, γ > 0 and the non-negativity of ∆, the gain l should satisfy l 1 ≤ l ≤ l 2 with l 1 , l 2 the positive roots of f (l), l 1,2 :=

12k 2 (k-1)λ± √ ∆ 2(2k-1) 2 λ 2 with ∆ = -16k 3 (k -1) 8λ 2 k 2 -(9λ 2 + 8λ 2 )k + 2λ 2 .
The proof can be concluded similarly to the last part of the proof of Theorem 1 in [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF].

Remark 1 We should note that there exist topologies and values of the gains k, l, that do not satisfy the conditions of our theorem but still ensure the consensus objective. Such scenario is, e.g the one considered in the simulations of [START_REF] Hong | Distributed observers design for leader-following control of multi-agent networks[END_REF]. This is intrinsically related to the selection of the Lyapunov function and, as such, it implies that another Lyapunov function can provide a better set of conditions.

2

  January 29, 2019 such that limt→∞ |x i (t) -x 0 (t)| = 0, lim t→∞ |v i (t) -v 0 (t)| = 0,namely, the agents can follow the leader (in the sense of both position and velocity).


Philippe Martinet

Keywords: servoing, vision-based control, robust control, automatic guided vehicles, active vision and sensor integration, visual tracking, and parallel architecture for visual servoing applications

Abstract-In this paper, it is shown that computer vision, used as an exteroceptive redundant metrology mean, simplifies the control of a Gough-Stewart parallel robot. Indeed, contrary to the usual methodology, where the robot is modeled independently from the control law which will be implemented, we take into account that vision will be used for control, from the early modeling stage. Hence, kinematic modeling and projective geometry are fused into a control-devoted projective kinematic model. Thus, a novel vision-based kinematic modeling of such a robot is proposed through the observation of its legs. Inspired by the geometry of lines, this model unifies and simplifies both identification and control. Indeed, it has a reduced parameter set, and allows us to propose a linear solution to its calibration. Using the same model, a visual servoing scheme is presented, where the attitudes of the nonrigidly linked legs are servoed, rather than the end-effector pose. Finally, theoretical results concerning the stability of this control law are provided.

Index Terms-Computer vision, line geometry, parallel kinematics, visual servoing.

I. INTRODUCTION

Parallel mechanisms are such that there exist several kinematic chains (or legs) between their base and their moving platform (also called the end-effector below). Therefore, they may exhibit a better repeatability [START_REF] Merlet | Parallel Robots[END_REF] than serial mechanisms, but not necessarily a better accuracy [START_REF] Wang | On the accuracy of a Stewart platform-Part I: The effect of manufacturing tolerances[END_REF], because of the large number of links and passive joints. There can be two ways to improve the accuracy. The first way is to perform a kinematic calibration of the mechanism, and the second one is to use a control law which is robust to calibration errors.

There exists a large amount of work on the control of parallel mechanisms (see [3] for a long list of references). In the focus of attention, Cartesian control is naturally achieved through the use of the inverse differential kinematic model (abusively called the robot Jacobian, since it can not be expressed as a matrix of partial derivatives of a map with respect to its coordinates), which transforms Cartesian velocities into joint velocities. It is noticeable that the inverse differential kinematic model of parallel mechanisms does not only depend on the joint configuration (as for serial mechanisms), but also on the end-effector pose.

Consequently, one needs to be able to estimate or measure the latter. As far as we know, all the effort has been put on the estimation of the end-effector pose through the forward kinematic model and the joint measurements. However, this yields much trouble, related to the fact that there is usually no analytic formulation of the forward kinematic model of a parallel mechanism. Hence, one numerically inverts the inverse kinematic model, which is algebraically defined for most of the parallel mechanisms. However, it is known [START_REF] Merlet | An algorithm for the forward kinematics of general 6 D[END_REF], [START_REF] Husty | An algorithm for solving the direct kinematics of general Gough-Stewart platforms[END_REF] that this numerical inversion requires high-order polynomial root determination, with several possible solutions (up to 40 real solutions [START_REF] Dietmaier | The Stewart-Gough platform of general geometry can have 40 real postures[END_REF] for a Gough-Stewart platform [START_REF] Gough | Universal tyre test machine[END_REF], [START_REF] Stewart | A platform with six degrees of freedom[END_REF]). Much of the work is thus devoted to solving this problem accurately and in real time (see, for instance, [START_REF] Zhao | Direct displacement analysis of parallel manipulators[END_REF]), or to designing parallel mechanisms with algebraic forward kinematic model [START_REF] Kim | Evaluation of a Cartesian parallel manipulator[END_REF], [START_REF] Gogu | Fully-isotropic T3Rl-type parallel manipulator[END_REF]. Alternately, one of the promising paths lies in the use of the so-called metrological redundancy [START_REF] Baron | The on-line direct kinematics of parallel manipulators under joint-sensor redundancy[END_REF], which simplifies the kinematic models by introducing additional sensors into the mechanism, and thus yields easier control [START_REF] Marquet | Contribution à l'étude de l'apport de la redondance en robotique parallele[END_REF].

Computer vision being an efficient way of estimating the endeffector pose [START_REF] Dementhon | Model-based object pose in 25 lines of code[END_REF], [START_REF] Lavest | Do we really need an accurate calibration pattern to achieve a reliable camera calibration[END_REF], it is a good alternative to use it for Cartesian control of parallel mechanisms. It can be done in three ways.

1) Vision as a Sensor:

The first one consists of computing the end-effector poses by vision, then in translating them into joint configurations through the inverse kinematic model, and finally servoing in the joint space. This scheme is rather easy to implement for serial mechanisms, provided that inverting the forward kinematic model can be done satisfactorily. The latter is straightforward for parallel mechanisms, since they usually have an algebraic inverse kinematic model. Similarly, one can consider computer vision as a contactless redundant sensor, as already stated in the context of parallel mechanism calibration [START_REF] Andreff | Vision-based kinematic calibration of an H4 parallel mechanism: Practical accuracies[END_REF], and use the simplified models based on the redundant metrology paradigm.

However, such schemes should be used carefully for parallel mechanisms, since joint control does not take into account the kinematic closures, and may, therefore, yield high internal forces [START_REF] Dasgupta | Force redundancy in parallel manipulators: Theoretical and practical issues[END_REF]. Moreover, there may exist several end-effector poses associated with a given joint configuration. Hence, a simple joint control may converge to the wrong end-effector pose, even though it converges to the correct joint configuration.

2) Visual Servoing: Second, vision can be additionally used to perform visual servoing [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]. Indeed, instead of measuring the end-effector pose and converting it into joint values, one could think of using this measure directly for control. Recall that there exist many visual servoing techniques, ranging from positionbased visual servoing (PBVS) [START_REF] Martinet | Vision based control law using 3D visual features[END_REF] (when the pose measurement is explicit) to image-based visual servoing (IBVS) [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] (when it is made implicit by using only image measurements). Most applications embed the vision system onto the end-effector to position the latter with respect to a rigid object whose accurate position is unknown, but one can also find applications with a fixed camera observing the end-effector [START_REF] Horaud | Visually guided object grasping[END_REF]. The interested reader is referred to [21] for a thorough and up-to-date state-of-the-art.

Visual servoing techniques are very effective, since they close the control loop over the vision sensor. This yields a high robustness to perturbations as well as to calibration errors. Indeed, these errors only appear in a Jacobian matrix, but not in the regulated error.

Essentially, visual servoing techniques generate a Cartesian desired velocity, which is converted into joint velocities by the robot inverse differential kinematic model. Hence, one can translate such techniques to parallel mechanisms, as in [START_REF] Koreichi | Visual servo control of a parallel manipulator for assembly tasks[END_REF]- [START_REF] Kallio | Three-dimensional position control of a parallel micromanipulator using visual servoing[END_REF] [for parallel robots with a reduced number of degrees of freedom (DOFs)], by observation of the robot end-effector and the use of standard kinematic models. It is rather easier than in the serial case, since the inverse differential kinematic model of a parallel mechanism is usually algebraic. Moreover, for parallel mechanisms, since the joint velocities are filtered through the inverse differential kinematic model, they are admissible, in the sense that they do not generate internal forces. More precisely, this is only true in the ideal case. However, if the estimated inverse differential kinematic model used for control is close enough to the actual one, the joint velocities will be closely admissible, in the sense that they do not generate high internal forces. The only difficulty for end-effector visual servoing of a parallel mechanism would come from the dependency of the inverse differential kinematic model to the Cartesian pose, which would need be estimated, but, as stated above, vision can also do that [START_REF] Dementhon | Model-based object pose in 25 lines of code[END_REF], [START_REF] Lavest | Do we really need an accurate calibration pattern to achieve a reliable camera calibration[END_REF]! Notice that this point pleads for PBVS rather than IBVS of parallel mechanisms, which is effectively the choice made in [START_REF] Koreichi | Visual servo control of a parallel manipulator for assembly tasks[END_REF]- [START_REF] Kallio | Three-dimensional position control of a parallel micromanipulator using visual servoing[END_REF].

From the above discussion, we thus highly recommend using visual servoing for parallel mechanism control.

3) A Novel Approach: However, the previous two ways consist solely of a simple adaptation of now classical control schemes, which, although probably very efficient, are not very innovative. Moreover, the direct application of visual servoing techniques assumes implicitly that the robot inverse differential kinematic model is given, and that it is calibrated. Therefore, modeling, identification, and control have small interaction with each other. Indeed, the model is usually defined for control using proprioceptive sensing only, and does not foresee the use of vision for control, then identification and control are defined later on with the constraints imposed by the model (Fig. 1). This is useful for modularity, but this might not be efficient in terms of accuracy, as well as in experimental setup time.

On the opposite side, a unified framework for modeling, identification, and control, apart from being definitely more satis- fying for the mind, would certainly open a path to higher efficiency. Indeed, instead of having identification and control being driven by the initial modeling stage, one could have a model taking into account the use of vision for control, and hence, for identification. To do so, it is necessary to fuse robot kinematics and projective geometry into a projective kinematic model (Fig. 2). Thus, we propose a novel third way to use vision, which gathers the advantages of redundant metrology and of visual servoing and avoids most of their drawbacks.

Moreover, observing the end-effector of a parallel mechanism by vision may be incompatible with its application. For instance, it is not wise to imagine observing the end-effector of a machining tool. On the opposite side, it should not be a problem to observe the legs of the mechanism, even in such extreme cases. Thereby one would turn vision from an exteroceptive sensor to a somewhat more proprioceptive sensor. This brings us back to the redundant metrology paradigm.

Parallel mechanisms are most often designed with slim and rectilinear legs. Thus, one is inclined to consider them as straight lines, as it was done for their kinematic analysis [START_REF] Merlet | Parallel Robots[END_REF], [START_REF] Husty | An algorithm for solving the direct kinematics of general Gough-Stewart platforms[END_REF] or kinematic calibration [START_REF] Andreff | Vision-based kinematic calibration of an H4 parallel mechanism: Practical accuracies[END_REF], [START_REF] Renaud | Kinematic calibration of parallel mechanisms: A novel approach using legs observation[END_REF]- [START_REF] Renaud | Combining end-effector and legs observation for kinematic calibration of parallel mechanisms[END_REF]. Therefore, the line geometry [START_REF] Plücker | On a new geometry of space[END_REF], [START_REF] Semple | Algebraic Projective Geometry[END_REF] is certainly the heart of the unification, all the more as line geometry is widely used in kinematic analysis [START_REF] Maccarthy | Introduction to Theoretical Kinematics[END_REF], [START_REF] Pottmann | Approximation in line space-Applications in robot kinematics and surface reconstruction[END_REF] and computer vision [START_REF] Faugeras | Three-Dimensional Computer Vision-A Geometric Viewpoint[END_REF], and has already been used for visual servoing [START_REF] Mahony | Visual servoing using linear features for under-actuated rigid-body dynamics[END_REF]- [START_REF] Mahony | Image-based visual servo control of aerial robotic systems using linear image features[END_REF].

Previous work on kinematic calibration [START_REF] Renaud | Kinematic calibration of parallel mechanisms: A novel approach using legs observation[END_REF]- [START_REF] Renaud | Combining end-effector and legs observation for kinematic calibration of parallel mechanisms[END_REF] already considered vision as a way to deliver contactless metrological redundancy. However, with the exception of [START_REF] Renaud | Combining end-effector and legs observation for kinematic calibration of parallel mechanisms[END_REF], the models that were calibrated remain classical. Indeed, vision was only used for sensing, and neither modeling nor control was questioned from the vision point of view. A first step in this direction was made in [START_REF] Andreff | Vision-based control of a Gough-Stewart parallel mechanism using legs observation[END_REF], where vision was used already at the modeling stage in order to derive a visual servoing scheme based on the observation of a Gough-Stewart parallel robot [START_REF] Gough | Universal tyre test machine[END_REF], [START_REF] Stewart | A platform with six degrees of freedom[END_REF]. However, full analysis of the model and its calibration was not addressed.

Consequently, the contribution of this paper is to present an original and unifying vision-based modeling and identification and control framework of parallel mechanisms by observing their legs with a camera fixed with respect to the base. This framework, inspired by line geometry, synthesizes into a single convergent framework previous work on identification [START_REF] Renaud | On vision-based kinematic calibration of n legs parallel mechanisms[END_REF] and on control [START_REF] Andreff | Vision-based control of a Gough-Stewart parallel mechanism using legs observation[END_REF], simplifying identification and going deeper into control analysis. It is introduced in the case of parallel mechanisms of the hexapod type, with illustration on a Gough-Stewart platform (Fig. 3).

The remainder of the paper is the following. Section II emphasizes the vision-based kinematics of the hexapod, which is the key point of the paper. A discussion of this model with respect to sensing follows in Section III. Then, Section IV recalls the differential geometry aspect of the leg observation and the control derived from it. Section V contains a closed-form linear solution to the identification of the model used for control. Finally, simulation and experimental results and the conclusion are given, respectively, in Sections VI and VII.

II. VISION-BASED KINEMATICS

A. Preliminaries

This paper assumes that the vision system is calibrated, which is not a strong hypothesis anymore, since accurate and easy-to-use camera calibration software can easily be found on the Web [38]. The classical approach for calibrating a camera is to use a calibration grid made of points [START_REF] Lavest | Do we really need an accurate calibration pattern to achieve a reliable camera calibration[END_REF], [START_REF] Brown | Close-range camera calibration[END_REF]- [START_REF] Faugeras | Camera calibration for 3D computer vision[END_REF]. However, in this case, where we aim at a unifying framework using line geometry, it may be more satisfactory for the mind to use the method proposed in [START_REF] Devernay | Straight lines have to be straight[END_REF]. Indeed, this method is based on the observation of a set of lines, without any assumption on their 3-D structure. As can be seen in Fig. 3, there are plenty of lines in the images observed by a camera placed in front of a Gough-Stewart parallel robot. Hence, the camera may be calibrated online without any experimental burden.

Since we plan to use line geometry as a central element for modeling, a representation for lines suited to control and identification is needed. Among the work on visual servoing from lines [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], [START_REF] Mahony | Visual servoing using linear features for under-actuated rigid-body dynamics[END_REF]- [START_REF] Mahony | Image-based visual servo control of aerial robotic systems using linear image features[END_REF], [START_REF] Malis | Intrinsics-free visual servoing with respect to straight lines[END_REF], [START_REF] Khadraoui | Vision based control in driving assistance of agricultural vehicles[END_REF], we prefer the so-called binormalized Plücker coordinates representation in [START_REF] Andreff | Visual servoing from lines[END_REF], which turns out to be coherent with kinematic modeling of parallel mechanisms.

In such a representation, a straight line in the oriented 3-D space [START_REF] Stolfi | Oriented Projective Geometry[END_REF] is modeled by the triplet where: • is the unit vector giving the spatial direction of the line; • is also a unit vector, and is a nonnegative scalar. They are defined by the cross-product

, where is any point on the line and depends on the reference point . Notice that using this notation, the well-known (normalized) Plücker coordinates [START_REF] Plücker | On a new geometry of space[END_REF], [START_REF] Pottmann | Approximation in line space-Applications in robot kinematics and surface reconstruction[END_REF] are the couple . An interesting property of this representation, concerning computer vision, is that if is chosen as the center of projection, then represents the image projection of the line, i.e., the equation of the image line verifies [START_REF] Merlet | Parallel Robots[END_REF] where and are the coordinates of a point in the image. The interpretation of the scalar is the orthogonal distance of the line to the center of projection.

Another interesting property of this representation is that for any point in the 3-D space, the dot product is the signed orthogonal distance of to the interpretation plane (i.e., the plane containing the projection center and the image line), defined by .

B. Kinematics of an Hexapod

Consider the hexapod in Fig. 3. It has six legs of varying length , attached to the base by spherical joints located in points , and to the moving platform (end-effector) by spherical joints located in points . 1 The inverse kinematic model of such an hexapod is given by [START_REF] Wang | On the accuracy of a Stewart platform-Part I: The effect of manufacturing tolerances[END_REF] expressing that is the length of vector . This model can be expressed in any Euclidean reference frame. Hence, it can be expressed in the base frame , in the end-effector frame , or in the camera frame . In the remainder and when needed, the reference frame used will be made explicit by a left superscript.

Let us consider , the binormalized Plücker coordinates of the line passing through and , oriented from to . Then, one trivially has

(3) (4)
From [START_REF] Merlet | Parallel Robots[END_REF], it is known that the inverse differential kinematic model of the hexapod, relating the end-effector Cartesian velocity (i.e., its kinematic twist) to the joint velocities is . . . . . . [START_REF] Husty | An algorithm for solving the direct kinematics of general Gough-Stewart platforms[END_REF] where is the center of the end-effector reference frame. Notice that the inverse differential kinematic model is written rather than , to clearly state that it has an algebraic expression, contrary to the inverse differential kinematic model of a serial mechanism.

C. Vision-Based Kinematics of a Hexapod

It has been noticed [START_REF] Merlet | Parallel Robots[END_REF] that the lines of the inverse differential kinematic model are the Plücker coordinates of the legs. Indeed, they are the wrenches [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems[END_REF] applied to the moving platform by the actuated legs. Under both interpretations, their expression depends on the chosen reference point. The advantage of taking this point on the mobile platform (for instance, ) is that in such a case, are constant, and only depends on

. Consequently, if one can measure or estimate (i.e., through a model) in the end-effector frame, one can easily convert , the end-effector Cartesian velocity expressed in the end-effector frame into joint velocities.

This measure can be done with a camera embedded onto the end-effector (i.e., in a first approximation

) and observing the legs (see Section III-B). In this case, the vision-based kinematics of the hexapod expressed in the end-effector frame are very simple

This formulation can apply to a classical visual servoing scheme with an embedded camera. Indeed, such a scheme generates, without loss of generality, a desired from the images of an externally fixed target. In practice, it may be awkward, since the camera should observe both the external target and all the legs. Alternately, several cameras could be used, but need be synchronized and calibrated with respect to each other.

In practice, it may be more convenient if the camera observing the legs is fixed to the base. Then, the reference frame associated to it is, without loss of generality, the base frame, and the kinematics of the hexapod do not express as simply as in the

where . Hence, it is necessary with this expression to estimate the end-effector orientation with respect to the base frame.

An alternate formulation is possible, which is somewhat less useful for standard Cartesian control. However, it is well suited to the observation of the legs only, and thereby to the control scheme proposed in Section IV. It consists of considering the mechanism in its dual operating mode: the end-effector is fixed, and the base moves with respect to it. Thus, we are interested in the inverse differential kinematic model relating the base Cartesian velocity expressed in the base frame to the joint velocities. By analogy with ( 6)-( 9

) ( 13
) 14
Notice the minus signs in ( 11) and ( 13), coming from the fact that in the permutation the direction of the legs has changed (Fig. 4). Notice also that now the inverse differential kinematic model is independent from the relative pose of the end-effector and base, since it only depends on the measured 's and the constant kinematic parameters 's. Notice also that ( 14) uses the fact that is the same for any point on the line.

D. Vision-Based Kinematics Expressed in the Camera Frame

In practice, the camera is not co-located with the base frame of the robot, nor share the same orientation in space. Consequently, the camera-to-base transformation differs from the identity, and the above model ( 11)-( 13) may be questioned.

Nevertheless, it is easy to show that one only needs to replace by in the super-and subscripts to convert the vision-based kinematics expressed in the base frame into vision-based kinematics expressed in the camera frame.

III. DISCUSSION

Here, the proposed model is discussed, with explicit regard to the sensing problem and to control. Indeed, no matter the control law, it will make use of some differential kinematic model, and thus will need to measure or estimate the 's. Moreover, the derivation of ( 11)-( 13) does not absolutely imply that computer vision should be used for that.

A. Why Vision Should Be Used Rather Than Joint Sensors

There are three manners to measure or estimate the 's. The first one, which, of course, we discard immediately, is to estimate the end-effector pose with respect to the base by numerical inversion of the inverse kinematic model, and then use [START_REF] Gogu | Fully-isotropic T3Rl-type parallel manipulator[END_REF] to obtain .

The second manner is to place joint sensors in the 's, so that they would deliver the pointing direction of the leg. This manner is valid, since the 's would be measured and not estimated through a delicate numerical procedure, as above. theless, is, in our opinion, still not the correct manner. Indeed, a joint sensor delivers a value expressed in its own reference frame. To convert this value in the base frame would require either an extremely accurate assembly procedure, or the identification of the relative orientation of each joint sensor frame with respect to the base frame. Moreover, additional joint offsets would need be identified.

Since the leading vector of a leg is essentially a Cartesian feature, we chose to estimate it by vision. Indeed, vision is an adequate tool for Cartesian sensing, and, following [START_REF] Renaud | On vision-based kinematic calibration of a Stewart-Gough platform[END_REF], if vision is also chosen for calibration, this does not add an extra calibration parameter. It will even be shown in Section V that using vision reduces the parameter set needed for control.

B. Cylindrical Leg Observation

Now the problem is to recover from the leg observation. It may be somehow tedious, although certainly feasible, in the case of an arbitrary shape. Hopefully, for mechanical reasons such as rigidity, most of the parallel mechanisms are not only designed with slim and rectilinear legs, but, even better, with cylindrical shapes.

Except in the degenerated case where the projection center lies on the cylinder axis, the visual edge of a cylinder is a straight line (Fig. 5). Consequently, it can be represented by its binormalized Plücker coordinates in the camera frame. Let us note and as the image projections of the two edges of a cylinder. These two vectors are oriented so that they point from the cylinder revolution axis outwards. Then, their expression is related to the binormalized Plücker coordinates of the cylinder axis (see Fig. 6) by ( 15) (16) where and , and is the cylinder radius.

From ( 15) and ( 16), it is easy to show that the leading vector of the cylinder axis, expressed in the camera frame, writes [START_REF] Dasgupta | Force redundancy in parallel manipulators: Theoretical and practical issues[END_REF] Notice that the geometric interpretation of this result is that is, up to a scale factor, the vanishing point of the two image edges, i.e., their intersection point in the image.

IV. VISION-BASED CONTROL

In this section, the control problem is addressed: from a given configuration of the hexapod legs observed by a camera attached to the base, how to reach a desired configuration?

Visual servoing is based on the so-called interaction matrix [START_REF] Chaumette | La commande des robots manipulateurs[END_REF] which relates the instantaneous relative motion between the camera and the scene, to the time derivative of the vector of all the visual primitives that are used through [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] where and are, respectively, the kinematic twist 2 of the camera and the scene, both expressed in . Then, one achieves exponential decay of an error between the current primitive vector and the desired one using a proportional linearizing and decoupling control scheme of the form [START_REF] Martinet | Vision based control law using 3D visual features[END_REF] where is used as a pseudocontrol variable. In this paper, we will similarly need to define a visual primitive, then form an error between its current value and its desired one, then relate in some way its time derivative to the actuation, and finally find a control relation between the error and the actuation.

A. Visual Primitive and Error

As foreseen above, the unit vectors will be used as visual primitives. Since these primitives are expressed in the 3-D space, the control scheme is close to a PBVS scheme. However, since the reconstruction step ( 17) is algebraic, it is, nevertheless, not far from

The visual primitives being unit vectors, it is theoretically more elegant to use the geodesic error rather than the standard vector difference. Consequently, the error grounding the proposed control law will be (20)

B. Interaction Matrix

Here, the time derivative of is related to the actuation. From (3), one immediately obtains (21) Inserting the interaction matrix associated with a 3-D point [START_REF] Martinet | Vision based control law using 3D visual features[END_REF] applied to the moving point [START_REF] Koreichi | Visual servo control of a parallel manipulator for assembly tasks[END_REF] where is the antisymmetric matrix associated with the cross product into (21), yields [START_REF] Kino | A motion control scheme in task oriented coordinates and its robustness for parallel wire driven systems[END_REF] It is interesting to see that both the base Cartesian velocity and the joint velocity vector appear in this expression, while also being linked to each other by the inverse differential kinematic model in [START_REF] Marquet | Contribution à l'étude de l'apport de la redondance en robotique parallele[END_REF]. This is certainly due to the existence of closed kinematic chains.

Nevertheless, using precisely the linking inverse differential kinematic model, one can exhibit a relationship between each and only. Indeed, each line of the inverse differential kinematic model in (13) rewrites as [START_REF] Kallio | Three-dimensional position control of a parallel micromanipulator using visual servoing[END_REF] Hence, one gets the following relationship:

(25) (26)
where is 3 6 and is obviously of rank 2.

Since , one can rewrite the above expression, using uniquely constant or easily measurable ( and ) quantities, as [START_REF] Renaud | Combining end-effector and legs observation for kinematic calibration of parallel mechanisms[END_REF] Necessary Condition 1: A minimum of three independent legs is necessary to control the end-effector pose, provided that there exists a diffeomorphism between the task space and the end-effector pose (i.e., provided that the observation of three legs allows determining the end-effector pose).

However, it should be borne in mind that this is solely a necessary condition, but that it does not hold any sufficient counterpart. The implication of this necessary condition is that the proposed approach might be used for parallel mechanisms similar to the Gough platform, but with less than six DOFs.

A global pseudointeraction matrix can then be obtained by stacking each individual interaction matrices . However, it is, in our opinion, an open question whether shall or shall not be considered as an interaction matrix. Indeed, in visual servoing, the various visual primitives are the image projections of objects in space that are rigidly linked to each other, while here, each of the legs is in relative motion with respect to the other ones. Nevertheless, effective control can be derived as shown in the following section.

C. Control Law

Let us choose a control such that decreases exponentially, i.e., such that [START_REF] Plücker | On a new geometry of space[END_REF] Then, introducing and , the combination of ( 20), [START_REF] Renaud | Kinematic calibration of parallel mechanisms: A novel approach using legs observation[END_REF], and (28) gives [START_REF] Semple | Algebraic Projective Geometry[END_REF] which yields [START_REF] Maccarthy | Introduction to Theoretical Kinematics[END_REF] The Cartesian pseudocontrol velocity is, hence [START_REF] Pottmann | Approximation in line space-Applications in robot kinematics and surface reconstruction[END_REF] and can be transformed into the control joint velocities using (12) [START_REF] Faugeras | Three-Dimensional Computer Vision-A Geometric Viewpoint[END_REF] where the hat means that only an estimate can be used.

Notice that since the joint values appear marginally in [START_REF] Renaud | Combining end-effector and legs observation for kinematic calibration of parallel mechanisms[END_REF], a control scheme which does not make any use of the joint values can be set up by using the median joint values [START_REF] Andreff | Visual servoing of a Gough-Stewart parallel robot without proprioceptive sensors[END_REF].

D. Stability

Let us consider the error function . Then, the proposed control rewrites as . Let us now consider the time derivative of under the latter control. It can easily be shown, from [START_REF] Semple | Algebraic Projective Geometry[END_REF] and inverting [START_REF] Baron | The on-line direct kinematics of parallel manipulators under joint-sensor redundancy[END_REF], that [START_REF] Mahony | Visual servoing using linear features for under-actuated rigid-body dynamics[END_REF] where encodes the actual motion of the end-effector under the effect of the joint actuation and all unmodeled phenomena.

As a consequence, an evident exponential stability condition of the error function is [START_REF] Andreff | Visual servoing from lines[END_REF] This condition is ensured if the following two conditions are satisfied:

In turn, these conditions are satisfied if the measured 's are sufficiently close to the real ones, and if the 's are accurately enough calibrated and if no singular configuration occurs. The singular configurations can be of two kinds: the much-studied kinematic singularities of the robot itself, or configurations such that is rank-deficient (that are left to future investigation).

V. VISION-BASED CALIBRATION

As seen above, the only kinematic parameters the control law depends on are the attachment points of the legs onto the base expressed in the camera frame and the joint offsets. The latter appear in two places in [START_REF] Renaud | Combining end-effector and legs observation for kinematic calibration of parallel mechanisms[END_REF]: under the form and as a gain. Considering the order of magnitude of and , one can neglect small errors on the joint offsets. Moreover, since the joints are prismatic, it is easy to measure their offsets manually with a millimetric accuracy. This is also highly sufficient to ensure that the gain is accurate enough. This means that, as far as control is concerned, there is no need for identifying the other usual kinematic parameters: attachment points onto the mobile platform and joint offsets.

In [START_REF] Renaud | On vision-based kinematic calibration of a Stewart-Gough platform[END_REF], a calibration procedure was proposed, using legs observation, where in a first step, the points were estimated in the camera frame, then in a second step, were expressed in the base frame, and finally, the others kinematic parameters where deduced. Essentially, only the first step of this procedure is needed.

This step is reformulated here in a more elegant way, using the binormalized Plücker coordinates of the cylinder edges [START_REF] Lavest | Do we really need an accurate calibration pattern to achieve a reliable camera calibration[END_REF], and minimizing an image-based criterion.

Assuming that the attachment point is lying on the revolution axis of the leg and referring again to Fig. 6, one obtains, for any leg and robot configuration [START_REF] Renaud | On vision-based kinematic calibration of n legs parallel mechanisms[END_REF] (38)

Stacking such relations for robot configurations yields the following linear system: . . . . . . [START_REF] Brown | Close-range camera calibration[END_REF] which has a unique least-square solution if there are at least two configurations with different leg directions. The calibration procedure is, hence, reduced to a strict minimum. To improve its numerical efficiency, one should only take care to use robot configurations with the larger angles between each leg direction. However, since the only appear in the pseudointeraction matrix, they do not require a very accurate estimation.

One could argue that, in practice, using a CAD model might be largely sufficient, except that a CAD model is not expressed in the camera frame, but in the base frame. Thus, one either has to estimate the base-to-camera transformation (which is a nonlinear problem), or to simply run our linear calibration scheme.

VI. RESULTS

The commercial DeltaLab Table de Stewart in Fig. 3

A. Image Noise

It is not immediate to model realistically the effect of image noise on detecting the visual edges of a cylinder. One can imagine to use the representation of an image line and to add noise to these two components. However, it is not certain that the two noise components are uncorrelated. Alternately, one can also determine the virtual intersections of the visual edge with the image border and to add noise to these intersections [START_REF] Bartoli | The 3D line motion matrix and alignment of line reconstructions[END_REF].

An alternative model is presented here, which takes advantage of the fact that an image line is essentially a unit vector. Thus, image noise will necessary result into a rotation of this unit vector. Consequently, one needs to define a random rotation matrix, that is to say, a rotation axis and a positive rotation angle. Preliminary tests showed that simply taking this axis as a uniformly distributed unit vector and this angle as a uniformly distributed positive scalar gives a realistic behavior.

This noise model has the advantages that it is easy to implement, it does not depend on the simulated image size, and is parametered by a single scalar (i.e., the maximal amplitude of the rotation angle).

To give an idea of how to choose this maximal amplitude, an error of about pixel on the extremities of a 300 pixel-long line segment yields a rotation angle of approximately 0.05 .

B. Control Validation

In all the simulations presented here and below, the initial configuration of the platform is the reference configuration where all the legs have minimal length. The goal configuration is obtained from this reference configuration by a translation by 10 cm along the axis of the platform (upward vertical) and a rotation of 15 around the axis, thus reaching the workspace limit.

The control law proposed in this paper was first tested without any noise for validation and analysis purposes.

In a first simulation, all the legs are used for control. Fig. 7 that the errors on each leg converge exponentially to zero, with a perfect decoupling in the task space. Moreover, the joint velocities (Fig. 8) have a smooth behaviour. Finally, Fig. 9 shows that the desired end-effector pose is reached and the trajectory of the legs in the image.

C. Calibration Validation

In order to quantify the calibration procedure, the simulated robot was calibrated by moving it into its 64 extremal configurations (i.e., each leg joint is extended to its lower limit, then to its upper limit). In each configuration, the visual edges of each leg were generated from the inverse kinematic model and added noise as described above, with a maximal amplitude of 0.05 . This calibration procedure was repeated 100 times. As a result, the median error on each attachment point is less than 1 mm. Recall that this result is obtained linearly without making any prior assumption on the mechanism geometry, nor on the camera-base transformation.

D. Realistic Simulation

Now, a more realistic simulation is presented, where the calibration is first performed using the 64 extremal configurations, then control is launched using the calibrated values. Noise is added as above during image detection, with amplitudes of 0.01 , 0.05 , and 0.1 .

It is noticeable that the calibration errors (in terms of maximal error on each of the component of each attachment point) is, respectively, of 0.5, 1.4, and 10 mm. Additionally, the median error of the convergence tails are, respectively, 0.1, 0.6, and 1.1 mm, while the maximal error are 0.6, 1.9, and 3 mm. Graphically (Fig. 10), the sum of the errors on each leg still decreases, yielding a potentially good robustness.

E. Experiments

The experimental robot has an analog joint position controller that we interfaced with Linux-RTAI. Joint velocity control is emulated through this position controller with an approximate 20 ms sampling period (this part is not yet running under RTAI, but only under standard Linux). Frame grabbing, line tracking, and numerical computation are performed using ViSP, an open C++ library for visual servoing [START_REF] Marchand | ViSP for visual servoing: A generic software platform with a wide class of robot control skills[END_REF]. The camera used was a 15 Hz, 1024 768, IEEE1394 camera with 4.8 mm focal length placed approximately 1 m from the base center (see Fig. 3). It also has to be noticed that the mechanism presents high frictious disturbances that have not yet been compensated for, since friction seems to depend nontrivially on the robot configuration. Hence, to overcome these disturbances, we implemented the visual servoing control with an adaptive gain, function of the controlled error norm: low at init, high near convergence.

The platform was asked to move from a configuration where all the legs are retracted to one where the two rear legs are stretched out (see Fig. 11). The resulting evolution of the controlled leg direction with respect to time is displayed in Fig. 12, showing the convergence of the control despite a high noise level.

VII. CONCLUSION

A novel approach was proposed for controlling a parallel mechanism using vision as a redundant sensor, adding a proprioceptive nature to the usual exteroceptive nature of vision. Within this approach, the standard tryptic "modeling, identification, and control" is reformulated, yielding a unified vision-based modeling, identification, and control framework. It is not only theoretically extremely elegant, but appears to own interesting properties. Indeed, stability conditions were exhibited, as well as po- tentially good robustness properties. It was validated and illustrated on a Gough-Stewart platform. Nevertheless, the results presented here should apply for any robot of the same kind with at least three length-actuated legs linked with passive revolute, universal, or spherical joints to the base and end-effector.

However, this paper is only the seed of a vast research domain. Indeed, there are several points to be addressed before a safe and satisfying implementation can be made on a real platform. First, this control does not take into account joint limit avoidance. This point is fundamental, since these limits can easily be reached and their avoidance may not be as trivial as for serial mechanisms. Second, this control assumes a detection of cylinder edges, which is known to be delicate in vision. Moreover, as can be seen in Fig. 11, using a single camera may yield visual self-occlusions of the hexapod. To avoid such occlusions, one may use several cameras, but this would require additional calibration, intercamera synchronization, and collaboration. Alternately, we plan to use a panoramic camera placed between the mechanism legs.

Third, since the control is essentially based on the direction of each leg, one may think of extracting the latter from the image of a generally shaped leg. Finally, we plan to extend this framework to any parallel mechanism, and we started to put effort into the development of high-speed vision tools to cope with the high velocities of parallel mechanisms [START_REF] Ait-Aider | Simultaneous pose and velocity measurement for high-speed robots[END_REF].

Fig. 1 .

 1 Fig. 1. Usual cascade from modeling to vision-based control.

Fig. 2 .

 2 Fig. 2. Simplified cascade from modeling to vision-based control using a projective kinematic model.

Fig. 3 .

 3 Fig. 3. Experimental set-up: A Gough-Stewart platform observed by a camera.

Fig. 4 .

 4 Fig. 4. Duality between the mobile end-effector mode and the fixed end-effector mode.

  ), i.e., by permutation of the roles of and , and of and (Fig.4), one obtains the vision-based kinematics of the hexapod expressed in the base frame

Fig. 5 .

 5 Fig. 5. Projection of a cylinder in the image.

Fig. 6 .

 6 Fig. 6. Construction of the visual edges of a cylinder. The cylinder is viewed from the top.

  was simulated and controlled experimentally. It is such that for all with 270 mm, 195 mm, , and the legs range are [345 mm, 485 mm].

Fig. 7 .

 7 Fig. 7. Errors (unit less) on each leg e e with respect to time (expressed as iteration number).

Fig. 8 .

 8 Fig. 8. Joint velocities (m/s) with respect to time.

Fig. 9 .

 9 Fig. 9. Top: Trajectory in space with initial (magenta, dashed) and desired (red, dash-dotted) position of the platform. Bottom: Image trajectories.

Fig. 10 .

 10 Fig. 10. Robustness to noise: sum of squares of the errors E E with a noise amplitude of 0 , 0.01 , 0.05 , and 0.1 .

Fig. 11 .

 11 Fig. 11. Composition of the desired (dashed platform in the CAD view, black in the camera image) and initial (solid platform in the CAD view, white in the camera image) configurations.

Fig. 12 .

 12 Fig. 12. Evolution of the controlled leg direction errors with respect to time.

Notice that the spherical joints imply an uncontrollable rotation of the legs around their axis, which has no effect on the end-effector position.

Note that there is a usual mistake in most of the papers dealing with visual servoing between screws and twists, refer, for instance, to[START_REF] Angeles | Fundamentals of Robotic Mechanical Systems[END_REF] and[START_REF] Gallardo | Dynamics of parallel manipulators by means of screw theory[END_REF].

This work was supported in part by the CPER Auvergne 2003-2005 program and in part by the European Commu-nity under the IP NEXT Project 011815. This paper was presented in part at the IEEE International