Gustavo Arechavaleta

Efraín López-Damian

José Luis Morales
email: jmorales@itam.mx

E López-

On the Use of Iterative LCP Solvers for Dry Frictional Contacts in Grasping

. The underlying idea of that method is to combine projected Gauss-Seidel iterations with subspace minimization steps. Gauss-Seidel iterations are aimed to obtain a high quality estimation of the active set. Subspace minimization steps focus on the accurate computation of the inactive components of the solution. Overall the new method is able to compute fast and accurate solutions of severely ill-conditioned LCPs.

We compare the performance of a modification of the iterative method of Morales et al with Lemke's algorithm on robotic object grasping problems.

I. INTRODUCTION

Contact dynamics is a well-known and important topic in physics-based animation and robotics among other disciplines. Video games, for example, are a good source of medium to large-scale contact problems. The common issue in these fields is the fast computation of contacts. Constraints, however, strongly characterize the application. In computer animation the constraints are mainly intended to produce visually plausible motions. Thus, simple friction models suffice to preserve realism in the animation. In robotics, say for problems regarding stable biped locomotion and grasping tasks, contact constraints must guarantee feasible robot motions in the presence of friction [START_REF] Miller | Implementation of multi-rigid-body dynamics within a robotic grasping simulator[END_REF].

It turns out that solving contact problems is a timeconsuming process mainly because at each time step of the simulation, performed by numerical integration, a linear complementarity problem has to be solved.

In this paper we address the problem of computing multirigid-body dynamics with frictional contacts. Specifically, we explore the use of the algorithm described in [START_REF] Morales | An algorithm for the fast solution of symmetric linear complementarity problems[END_REF] to determine the contact forces in the presence of friction between colliding rigid bodies.

In addition to the frictionless contact problem, in a wide variety of applications it is desirable to consider a friction model (e.g. the Coulomb's friction model). This fact increases the difficulty of the problem since there is coupling between the friction (tangential) and normal force components. In practice, the most common and widely used model of friction is the linearized version of the Coulomb friction law by a polyhedral approximation of the friction cone at each contact point [START_REF] Anitescu | Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementary Problems, ser. Nonlinear Dynamics[END_REF], [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and coulomb friction[END_REF]. Nevertheless, there exist other ways for approximating the Coulomb friction model [START_REF] Lacoursiere | Splitting methods for dry frictional contact problems in rigid multibody systems: preliminary performance results[END_REF].

A consistent and solvable mathematical model of dry frictional contacts has been proposed in [START_REF] Anitescu | Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementary Problems, ser. Nonlinear Dynamics[END_REF]; for variants based on LCP formulations see for example [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and coulomb friction[END_REF], [START_REF] Anitescu | Time-stepping for threedimensional rigid body dynamics[END_REF], [START_REF] Al-Fahed | Hard and soft fingered robot grippers, the linear complementary approach[END_REF]. Alternative strategies include spring and damper systems [START_REF] Mirtich | Impulse based simulation of rigid bodies[END_REF], although they tend to be unstable.

Among the wide range of LCP contact formulations either acceleration based models [START_REF] Trinkle | On dynamic multi-rigid-body contact problems with coulomb friction[END_REF], [START_REF] Pang | Complementary formulations and existence of solutions of dynamics multi-rigid-body contact problems with coulomb friction[END_REF], [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF], [START_REF] Baraff | Issues in computing contact forces for non-penetrating rigid bodies[END_REF], or velocity stepping models [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and coulomb friction[END_REF], [START_REF] Anitescu | Time-stepping for threedimensional rigid body dynamics[END_REF], [START_REF] Anitescu | Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementary Problems, ser. Nonlinear Dynamics[END_REF], the most common method to solve the LCP is still that of Lemke [START_REF] Lemke | Bimatrix equilibrium points and mathematical programming[END_REF], [START_REF] Sargent | An efficient implementation of the lemke algorithm and its extension to deal with upper and lower bounds[END_REF]. Pivotal methods are efficient for small problems but they do not exploit sparsity effectively and may suffer from numerical instability. Therefore they are not suitable for large-scale applications, or when the LCPs are ill-conditioned; see for example [START_REF] Morales | Computational experience with several methods for large-scale convex quadratic programming[END_REF]. In the robotics literature, some improvements on Lemke's algorithm have been proposed in [START_REF] Lloyd | Fast implementation of lemke's algorithm for rigid body contact simulation[END_REF], [START_REF] Yamane | A numerically robust lcp solver for simulating articulated rigid bodies in contact[END_REF] to overcome some of its inherent deficiencies.

Iterative LCP solvers, including the projected Gauss-Seidel (PGS) algorithm, have been utilized to compute frictional contacts problems (e.g. [START_REF] Catto | Iterative dynamics with temporal coherence[END_REF], [START_REF] Jourdan | A gauss-seidel like algorithm to solve frictional contact problems[END_REF], [START_REF] Liu | Computation of three-dimensional rigid-body dynamics with multiple unilateral contacts using time-stepping and gauss-seidel methods[END_REF]) since they are numerically robust. However, they can be very slow for medium to large-scale and ill-conditioned problems. In grasping tasks, the contact LCPs are ill-conditioned in most of the cases. These evidences motivated us to use a modification of the algorithm described in [START_REF] Morales | An algorithm for the fast solution of symmetric linear complementarity problems[END_REF]. In the forthcoming discussion we will refer to this method as PGS-SM.

This paper is organized as follows. In Section II we recall the general form of a SLCP. Section III gives a short description of the two stages of PGS-SM. In Section IV we review the standard mixed LCP formulation for multi-rigidbody dynamics with dry frictional contacts. We also describe the approximation of the friction model used to formulate the SLCP contact problem. Section V presents a series of numerical experiments and simulations that illustrate the effectiveness of PGS-SM. Finally, we provide in Section VI some concluding remarks.

II. STRUCTURE OF THE SLCP

The form of a SLCP that will be used in the paper is:

0 ≤ u ⊥ Au + a ≥ 0 (1)
where the unknown of the problem is u, the matrix A and the vectors a are given. We assume that A is n × n symmetric positive definite.

We recall that (1) corresponds to the first-order optimality conditions of the bound constrained quadratic program

min u φ(u) = 1 2 u T Au + a T u s.t. u ≥ 0. (2)
The optimization literature on methods for solving bound constrained problems is wide. Gradient projection methods and interior-point methods are among the best techniques to solve this kind of problems. However, the computational experience reported in [START_REF] Morales | An algorithm for the fast solution of symmetric linear complementarity problems[END_REF] indicates that if speed is a hard constraint then PGS-SM is a strong competitor.

III. THE ITERATIVE ALGORITHM

The underlying idea of the iterative algorithm is to alternate two strategies. First, a few iterations of the projected Gauss-Seidel (PGS) method are applied to identify an estimate of the optimal active set. Then, the algorithm performs consecutive subspace minimization (SM) steps to improve this estimate and to compute accurate values of the variables in the inactive set. The algorithm returns to the PGS phase to check the convergence tests, if satisfied the algorithm exits with the solution; otherwise it continues to the SM phase.

A. Projected Gauss-Seidel Method

We briefly outline one iteration of the PGS iterative method. The k + 1-th approximation can be defined in scalar form as follows Algorithm Given the k-th approximation u k ≥ 0.

for i = 1, 2, ..., n

∆u i = A -1 ii (a i + j<i A ij u k+1 j ) u k+1 i = max(0, u k i -∆u i , ). (3)
end Observe that, since A is symmetric positive definite, the iteration is well defined.

B. The PGS-SM Method

The algorithm can be formally described as follows.

Algorithm PGS-SM. Choose k PGS , k SM the maximum number of PGS and subspace minimization iterations respectively; choose tol > 0; choose an initial approximation u 0 ≥ 0.

Begin • PGS-stage. Perform k PGS iterations of PGS. if ||u k+1 -u k || ∞ ≤ tol then exit otherwise -Set u ← u k+1
-Form the approximation to the active set

A = { i | u i = 0 } -Form the corresponding reduced problem Âû + â = 0, û ≥ 0 (4)
by eliminating from the original problem (1) rows and columns whose indices belong to A.

• SM-stage.

Set k ← 0. while k < k SM 1) Solve the linear system Âû + â = 0 with any method suitable for symmetric positive matrices. 2) if û ≥ 0 then return to PGS-stage. otherwise -Compute the projection û ← max(0, û).

-Form a new reduced system Âû + â = 0.

-

Set ← k + 1. end Return to PGS-stage End
Observe that, due to the equivalence between the SLCP and the bound constrained optimization problem, all reduced problems are well defined.

IV. THE DRY FRICTIONAL CONTACT SLCP

We apply the multi-rigid-body contact formulation described in [START_REF] Anitescu | Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementary Problems, ser. Nonlinear Dynamics[END_REF], [START_REF] Miller | Implementation of multi-rigid-body dynamics within a robotic grasping simulator[END_REF], whose modeling strategy allows both unilateral (contacts) and bilateral (joints) constraints between rigid bodies. This is a time-stepping scheme, based on velocity, that makes use of an approximation of the Coulomb friction law.

A. Dynamic contact Model

We start from the Newton equation of motion

M v = F c + F e , (5)
where M is a block-diagonal symmetric positive definite matrix that collects the mass and rotational inertia of each body; v represents the acceleration;

F c = f j + f l + f n + f t
is the sum of joint, joint limit, normal and friction constraint forces. F e represents both Coriolis and external forces.

To solve the multi-body dynamics it is common to use an explicit Euler integration method over a time step from t to t + 1 for the equation of motion

M (v t+1 -v t ) = h(F c + F e ), ( 6
)
where h is the integration step length. We consider pairwise constraints between rigid bodies. The non-penetration constraints and frictional forces are expressed as inequalities that can be used to formulate the problem with complementarity conditions. These conditions are defined as

0 ≤ N T v t+1 ⊥ ĉn ≥ 0 (7) 0 ≤ DT v t+1 + E λ ⊥ β ≥ 0 (8) 0 ≤ μ ĉn -E T β ⊥ λ ≥ 0 (9)
where ĉn = [c

n ... c

(n)

n ] T and β = [β (1) ... β (nd) ] T are the normal and friction force components; n is the number of contacts; d is the number of friction directions at contact j. N = [N (1) ... N (n) ] and D = [D (1) ... D (nd) ] are the constraint matrices for normal and friction forces. Each column of N denotes the contact normal; d columns of D span the tangent space of friction forces at contact j; λ = [λ (1) ... λ (n) ] T is a vector of Lagrange multipliers; μ = diag[µ (1) ... µ (n) ] is a diagonal matrix of friction coefficients; E is a matrix defined as

E = diag[1], 1 = [1 ... 1] T ∈ R d .
The condition [START_REF] Al-Fahed | Hard and soft fingered robot grippers, the linear complementary approach[END_REF] serves to avoid interpenetration. The equations ( 8) and ( 9) are the complementary conditions of the linearized cone of Coulomb friction law ||f t || ≤ µ||f n ||.

B. Joint Constraints

As we will show in Section V, we have integrated the PGS-SM algorithm into the open source package, Graspit [START_REF] Miller | Graspit!: A versatile simulator for grasp analysis[END_REF]. The current representation of joints between two connected rigid bodies has the form of equality constraints on their relative velocities. Although, given that we do not constrain the relative position of the bodies, the joint connections occasionally break over time due to small drifts in the system. This is corrected by applying

ĴT v t+1 = j j = -ζ ∆P h , (10)
where Ĵ = [J (1) ... J (m) ]; m is the number of joints; ∆P represents the positional errors between connected bodies; ζ takes values in [0, 1] for controlling the error corrections. We then define the constraint impulses as follows

hf j = Ĵ ĉj , (11)
where ĉj = [c

... c

(m) j

] T is the vector of magnitudes of joint constraint impulses. In addition, we should consider that joints have mechanical limits, thus we can formulate another inequality constraint similar to the contact normal constraint as:

LT v t+1 ≥ 0 ⊥ ĉl ≥ 0, (12)
where L = [L (1) ... L (q) ], ĉl = [c

l ] T are the required impulses to keep the joints within their bounds and q is the number of reached joint limits.

C. Constructing the SLCP

Considering the constraints ( 7)-( 12) mentioned above and adding the equation of motion [START_REF] Anitescu | Time-stepping for threedimensional rigid body dynamics[END_REF] to the complementarity conditions, we obtain the following mixed LCP [START_REF] Anitescu | Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementary Problems, ser. Nonlinear Dynamics[END_REF]:

  M -Ĵ -H ĴT 0 0 H T 0 N     v t+1 ĉj z   +   -a -j b   =   0 0 w   (13) 0 ≤ z ⊥ w ≥ 0, where H = [ L N D 0], z = [ ĉl ĉn β λ] T ; a = M v t + hk; b = [0 0 0 0] T ; w = [ε ρ σ γ] T ; N has the form N =   0 0 0 0 0 E μ -E T 0   . ( 14
)
Solving for v t+1 and ĉj in [START_REF] Lemke | Bimatrix equilibrium points and mathematical programming[END_REF] we obtain a LCP of the form

(G + N )z + g = w, G = H T M -1 H- H T M -1 Ĵ( ĴT M -1 Ĵ) -1 ĴT M -1 H g = b + H T M -1 a- H T M -1 Ĵ( ĴT M -1 Ĵ) -1 ĴT M -1 a. (15)
The LCP can be solved using Lemke pivot-based algorithm [START_REF] Lemke | Bimatrix equilibrium points and mathematical programming[END_REF]- [START_REF] Yamane | A numerically robust lcp solver for simulating articulated rigid bodies in contact[END_REF], [START_REF] Cottle | Complementary pivot theory of mathematical programming[END_REF]. However, using this formulation, it is not possible to directly apply the PGS-SM algorithm since the matrix G + N is not positive definite. Only when N = 0, G is positive semi-definite and symmetric. To force the positive definiteness of G, and without losing the physical meaning of the formulation, we add to it a square diagonal matrix CF M (called constraint force mixing). The positive diagonal elements of CF M allow soft constraints.

To avoid N = 0, we replace the friction constraints ( 8) and ( 9) by the following complementarity conditions [START_REF] Lacoursiere | Splitting methods for dry frictional contact problems in rigid multibody systems: preliminary performance results[END_REF]:

Dv -ṽ+ + ṽ-= 0 0 ≤ β -β-⊥ ṽ+ ≥ 0, 0 ≤ β+ -β ⊥ ṽ-≥ 0 ( 16
)
where ṽ+ and ṽare the positive and negative components of the tangential velocity; β+ , βare the lower and upper bounds of β respectively. In this case, the matrix D is composed by two perpendicular directions, D j1 , D j2 . Thus, D is on the tangent plane at contact j and perpendicular to f n . With these parameters, we can define the size of a friction box by imposing a fixed bound on the tangential forces -µ cn ≤ β ji ≤ µ cn for i = 1, 2, cn is an estimate of c n (e.g. the last LCP solution before the current integration step). According to the above considerations, the system still has the form [START_REF] Lemke | Bimatrix equilibrium points and mathematical programming[END_REF]

but now H = [ L N D], z = [ ĉl ĉn β] T , a = M v t + hk, b = [0 0 0] T , w = [ε ρ σ] T
and N is the null matrix. Finally, we solve for v t+1 and ĉj to get [START_REF] Morales | Computational experience with several methods for large-scale convex quadratic programming[END_REF] with the complementarity conditions in [START_REF] Lloyd | Fast implementation of lemke's algorithm for rigid body contact simulation[END_REF] such that:

0 ≤ ĉl ⊥ ε ≥ 0, 0 ≤ ĉn ⊥ ρ ≥ 0, 0 ≤ β -β -⊥ σ ≥ 0, 0 ≤ β + -β ⊥ σ ≥ 0. (17)
Finding the LCP solution vector z, we then obtain

v t+1 = M -1 ( Ĵ ĉj + Hz + a) ĉj = -( ĴT M -1 Ĵ) -1 ĴM -1 (Hz + a)+ ( ĴT M -1 Ĵ) -1 j . V. SIMULATION EXPERIMENTS
In our experiments we compare three LCP methods to solve robotic object grasping with frictional contacts: the standard PGS method (III-A), the adapted PGS-SM method (III-B), and the standard Lemke method provided by Graspit [START_REF] Miller | Graspit!: A versatile simulator for grasp analysis[END_REF].

The experiments have been performed on a AMD Turion 64 X2 dual-core 3.0 GHz. Graspit was compiled using GNU C++. It includes a collision detection library called PQP [START_REF] Group | Pqp -a proximity query package[END_REF] and performs basic linear algebra operations with BLAS [START_REF] Blackford | An updated set of basic linear algebra subprograms (blas)[END_REF]. The PGS and PGS-SM algorithms have been implemented in Fortran 77 package using double precision and sparse linear algebra. The subspace minimization is carried out with preconditioned GMRES for nonsymmetric linear systems [START_REF] Saad | Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. Both algorithms were compiled using GNU f77.

The PGS and PGS-SM methods have been adapted for frictional contact problems by adding lower and upper bounds on the variables. We have integrated both methods in Graspit and calling them via a standard Fortran to C library.

A. Grasping Scenarios

The implementations of PGS and PGS-SM were tested using some scenarios provided by Graspit.

In the 3D grasping simulator, we defined an initial configuration of a robotic hand that can be attached to a Puma arm or not, and with the fingers opened. We also defined the initial configuration of a movable object on a static table and close to the robotic hand (see Fig. 1.(a), Fig. 2.(a)). We can see that objects are already in contact with the table surface, generating friction cones at these contacts.

We then activated the multi-rigid-body dynamics process considering both the robot and the object. At this stage, the hand started to close their fingers. After some iterations of the numerical integration process the fingers collide with the object producing some contacts between any part of the hand (fingers' phalanges and/or the palm), the object and the table. Depending on the pre-shape of the hand, the fingers do not always touch the object at the same time, causing the object to slip unless the grasp can balance the forces to finally carrying the object to a stable pose, see Figs. 1.(b)-(d)). The simulation is stopped by the user.

In object manipulation and more precisely, in grasp analysis and planning, a fundamental problem is the computation of stable grasps. This allows the robot to perform a series of tasks (e.g. pick and place) in human environments. Most of the planners only solve the problem from a geometric and kinematic point of view giving good results. However, it is necessary to incorporate dynamics in the analysis of grasping for practical applications in robotics. One of the most important elements of this analysis is the frictional contact model between the object and hand.

A grasp could seem to be stable but when dynamic effects are considered, we notice that it is not. We show in Fig. 2 four instances of a grasping task. We have a plastic hand while the object is made of glass, thus, the coefficient of friction between these materials is µ = 0.2. We observe that contact points appear with their respective frictional cones when the fingers touch the object (see Fig. 2.(b)). These contact points change in time due to the applied forces of the hand (see Fig. 2.(c)). In this example, the solutions of the multi-rigid-body dynamics in contact have been solved by the PGS-SM algorithm at every iteration. At the end of the simulation, we notice that even if the initial grasp seems to be good and the finger velocities are not so large, the hand cannot grasp the object (see Fig. 2(d)). This behavior frequently occurs for small coefficients of friction allowing slippage, which is a logic physical consequence. We perform the same experiment but changing only the hand material to rubber. In this case, the coefficient of friction is µ = 1.0. We see in Fig. 3.(a) the frictional cones with µ = 0.2 corresponding to wood table and glass cup while Fig. 3.(b) illustrates when the frictional cones between the hand (rubber) and the object. After some iterations the robot is capable to grasp the object stably.

B. Comparing Performances

To illustrate the behavior of the PGS-SM algorithm, we have conducted several comparisons between Lemke, PGS and PGS-SM in terms of computation time and number of contacts (size of the problem). We tested different examples of grasping tasks and we executed 3 times the same task invoking different LCP solver each trial.

At each time step throughout the simulation, we measured the CPU time required to find a solution of the frictional contact LCP.

In these experiments, before calling the PGS-SM method, we should specify the input parameters k PGS and k SM . After some tunning phase, we identified that k PGS ∈ [START_REF] Anitescu | Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementary Problems, ser. Nonlinear Dynamics[END_REF][START_REF] Anitescu | Time-stepping for threedimensional rigid body dynamics[END_REF] and k SM ∈ [START_REF] Morales | An algorithm for the fast solution of symmetric linear complementarity problems[END_REF][START_REF] Anitescu | Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementary Problems, ser. Nonlinear Dynamics[END_REF] are the adequate ranges. We set the stopping tolerance to 1.0e -8 and the matrix CF M = diag [1.0e -11]. The reason for these values are only for the numerical precision problem in computing operations when exact zero value is considered, this does not modify the algorithms or affect their performance.

We observed that the performance of the Lemke method was inefficient with respect to PGS and PGS-SM methods when the number of contacts increase. Note that the size of We see that the final grasp is stable unlike the previous case.

the problems depends on the number of contacts and joint limits reached. From 1 to 3 contacts, there are not significant difference between the 3 methods. From 4 to 12 contacts, the CPU time varied from 79 to 11820 microseconds. Table I shows the minimum, the maximum and the mean CPU time for Lemke, PGS and PGS-SM respectively. These quantities show the advantage of the adapted PGS-SM. We noted that the PGS method was able to identify the active set in few iterations as it is pointed out in [START_REF] Morales | An algorithm for the fast solution of symmetric linear complementarity problems[END_REF]. The adapted PGS-SM method is quite efficient due to the right balance between identifying the active set (PGS iterations) and moving toward optimality (subspace minimization cycles). In our experiments the PGS-SM algorithm only did from 3 to 6 PGS iterations and used one subspace minimization cycle to get the solution. The PGS-SM is also significantly faster than Lemke in most cases.

The frictional contact problems in our experiments can be classified as small but ill-conditioned problems. Because the size of the problems was small, the performance of Lemke algorithm was acceptable. The PGS-SM method, on the other hand, continue to be efficient for solving large-scale and ill-conditioned LCPs. It has been tested with large-scale frictionless contact problems extracted from the open source package Open Dynamic Engine (see [START_REF] Morales | An algorithm for the fast solution of symmetric linear complementarity problems[END_REF] for details).

These preliminary results indicate that robotic frictional contact problems would be gained computation time by using PGS-SM method without regarding the size of illconditioned LCPs.

Simulation sequences of the examples presented in this work can be found at: http://www.cinvestav.edu.mx/saltillo/robotica/garechav/grasping/

VI. CONCLUSIONS

In this paper we used an adaptation of the two-stage iterative algorithm proposed in [START_REF] Morales | An algorithm for the fast solution of symmetric linear complementarity problems[END_REF] to deal with ill-conditioned dry frictional contact problems in grasping tasks. We formulate contact dynamics in a similar manner as [START_REF] Anitescu | Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementary Problems, ser. Nonlinear Dynamics[END_REF].

The PGS-SM method allows us to solve constrained motions for robotic hands that have complex dynamics by applying appropriate iterative methods. Bounds in frictional variables are utilized to guarantee convergence to the grasping problems. We compute the frictional bounds by using the information of the previous iteration. As a result, the grasp motion has been successfully verified by simulations. We tested our approach with different friction coefficients and we compared the resulting behaviors with the solutions computed using Lemke algorithm. We are still working to find different ways to estimate the bounds. Nevertheless, our preliminarily results show that the PGS-SM method can be successfully applied in practice.

Other challenging problems arising in Humanoid robotics and motion planning for several and highly articulated mechanisms can be viewed as complex contact problems. For instance, the problem of computing whole-body stable object manipulation [START_REF] Hsiao | Imitation learning of whole-body grasps[END_REF], [START_REF] Yoshida | Whole-body motion planning for pivoting based manipulation by humanoids[END_REF]. Our current work is focusing on validating the use of the PGS-SM method in these classes of robotic problems.

  Fig.1. The figures show the results for one of the experiments. We used a five-fingered metal hand and a rubber ball object, after some iterations the robot picked the object up stably.

Fig. 2 .

 2 Fig.2. The figure shows a three-fingered Barrett hand (plastic) and a glass object with µ = 0.2. The fingers make contact with the object applying forces at the contact points. Due to the small friction between materials, the grasp cannot resist and the object slips.

Fig. 3 .

 3 Fig.3. The figure shows the experiment with the same three-fingered hand but now it is made of rubber and the object is made of glass, with µ = 1.0. We see that the final grasp is stable unlike the previous case.

TABLE I CPU

 I TIME COMPLEXITY (IN MICROSECONDS).

		Min	Max	Mean
	3-fingered hand
	-Lemke	146	7050	672
	-PGS	97	1332	207
	-PGS-SM	84	355	145
	5-fingered hand
	-Lemke	220 11820 1885
	-PGS	95	3348	238
	-PGS-SM	79	452	155

VII. ACKNOWLEDGMENTS

The first two authors acknowledge the financial support of the National Council of Science and Technology (CONACyT) within the project No. 84855. J.L. Morales was supported by Asociación Mexicana de Cultura A.C.


Remco I Leine
email: remco.leine@imes.mavt.ethz.ch

Ueli Aeberhard
email: ueli.aeberhard@imes.mavt.ethz.ch

A WEAK FORM OF HAMILTON'S PRINCIPLE AS VARIATIONAL INEQUALITY

Keywords: unilateral constraint, variational inequality, variational principle, principle of d'Alembert-Lagrange

The classical form of Hamilton's principle holds for conservative systems with perfect bilateral constraints. In this paper we derive Hamilton's principle for perfect unilateral constraints (involving impulsive motion) using weak variations. The resulting principle has the form of a variational inequality in Hilbert space.

Introduction

Classical Analytical Mechanics is concerned with mechanical systems with perfect bilateral (mostly holonomic) constraints [START_REF] Papastavridis | Analytical Mechanics: A Comprehensive Treatise on the Dynamics of Constrained Systems: for Engineers[END_REF] and is closely related with the calculus of variations as many principles of statics and dynamics are formulated in terms of variational problems, e.g. the principle of virtual work and Hamilton's principle (see for instance [START_REF] Lanczos | The Variational Principles of Mechanics[END_REF]). Unilateral constraints, which are basically inequality constraints, are completely ignored in classical Analytical Mechanics because inequalities are not discussed by the classical calculus of variations and, also, because unilateral constraints in dynamics lead to shocks with discontinuities in the velocity. The mathematical tools to handle unilateral problems in statics, and later also in dynamics, have only been developed since the last four decades. The field of Nonsmooth Dynamics is now rapidly developing. We refer the reader to the textbooks [START_REF] Leine | Dynamics and Bifurcations of Non-Smooth Mechanical Systems[END_REF][START_REF] Glocker | Set-Valued Force Laws, Dynamics of Non-Smooth Systems[END_REF][START_REF] Brogliato | Nonsmooth Mechanics[END_REF].

Variational problems involving convex inequality constraints are described by variational inequalities and where first introduced by Hartman and Stampacchia [START_REF] Hartman | On some nonlinear elliptic differential equations[END_REF] to study partial differential equations. The applicabil-ity of the theory has since been expanded to include problems from mechanics, finance, optimisation and game theory. References on variational inequalities can be found in the standard reference [START_REF] Kinderlehrer | An Introduction to Variational Inequalities and their Applications[END_REF] or in the more recent books (Goeleven et al., 2003a;Goeleven et al., 2003b). Variational problems with non-convex inequality constraints are described by hemi-variational inequalities and are discussed in the seminal work of Panagiotopoulos [START_REF] Panagiotopoulos | Hemivariational Inequalities: Applications in Mechanics and Engineering[END_REF]. Moreover, (hemi-)variational inequalities are related to (non-)convex superpotentials through the subderivative known from Nonsmooth Analysis, see the work of Moreau [START_REF] Moreau | La notion de superpotentiel et les liaisons unilatérales en élastostatique[END_REF]) and Panagiotopoulos [START_REF] Panagiotopoulos | Non-convex superpotentials in the sense of F.H. Clarke and applications[END_REF].

The extension of classical Analytical Mechanics to perfect unilateral constraints asks for a reformulation of the variational principles of mechanics in terms of variational inequalities. The principle of d'Alembert-Lagrange in inequality form has been discussed in [START_REF] Panagiotopoulos | Analytical mechanics: Addendum I. Inequality constraints with elastic impacts. The convex case[END_REF][START_REF] Panagiotopoulos | Inequality constraints with elastic impacts in deformable bodies. The convex case[END_REF][START_REF] Goeleven | Inequality forms of d'alembert's principle in mechanics of systems with holonomic unilateral constraints[END_REF][START_REF] Goeleven | Dynamic hemivariational inequalities and their applications[END_REF]May, 1984a;[START_REF] May | Hamilton's principle as substationarity principle[END_REF]. Various forms of the principle of Hamilton as variational inequality form can be derived from the principle of d'Alembert-Lagrange as variational inequality [START_REF] Panagiotopoulos | Analytical mechanics: Addendum I. Inequality constraints with elastic impacts. The convex case[END_REF][START_REF] Panagiotopoulos | Inequality constraints with elastic impacts in deformable bodies. The convex case[END_REF].

The current paper puts one existing form of the principle of Hamilton as variational inequality, which has already been derived in [START_REF] Panagiotopoulos | Analytical mechanics: Addendum I. Inequality constraints with elastic impacts. The convex case[END_REF][START_REF] Panagiotopoulos | Inequality constraints with elastic impacts in deformable bodies. The convex case[END_REF], within the context of weak and strong extrema [START_REF] Cesar | Necessary conditions and sufficient conditions of weak minimum for solutions with corner points[END_REF][START_REF] Troutman | Variational Calculus and Optimal Control: Optimization with Elementary Convexity[END_REF]. It therefore becomes much more clear how these different forms of Hamilton's principle have to be understood. Furthermore, the weak form does not impose any requirements on the energy dissipation (or conservation).

2 Principle of d'Alembert-Lagrange as variational inequality In this section we discuss the principle of d'Alembert-Lagrange in inequality form [START_REF] Goeleven | Inequality forms of d'alembert's principle in mechanics of systems with holonomic unilateral constraints[END_REF][START_REF] Panagiotopoulos | Analytical mechanics: Addendum I. Inequality constraints with elastic impacts. The convex case[END_REF]. Consider a pointmass m with position r ∈ R 3 which is subjected to a constraint force R and an external force F . The principle of virtual work states that the virtual work vanishes for all virtual displacements δr, i.e. δW = (mr -F -R) T δr = 0 ∀δr.

(1)

The virtual displacements δr = r *r are infinitely small differences of arbitrary virtual positions r * and the actual position r, such that time is kept constant. We first consider the case where the constraint is a bilateral holonomic scleronomic constraint g(r) = 0 with g ∈ C 1 (R 3 , R), in other words the position of the point-mass is constrained to the manifold K = {r ∈ R 3 | g(r) = 0} (see Figure 1). Virtual displacements δr are admissible with respect to the constraint if they belong to the tangent space T K (r) = {z ∈ R 3 | ∂g ∂r z = 0}. If the virtual work of the constraint force vanishes for all virtual displacements which are admissible with respect to the constraint, i.e.

R T δr = 0 ∀δr ∈ T K (r), (2)
then we speak of a perfect bilateral constraint. A perfect constraint force is therefore normal to the constraint manifold in the sense that

-R ∈ N K (r), (3)
where N K (r) is the set of all vectors which are normal to the tangent space T K (r). Equations ( 2) and (3) are equivalent if K is a manifold. Instead of taking (2) as definition of a perfect constraint, we therefore could also have chosen to take the normality condition (2) of the constraint force as definition. The assumption of perfect bilateral constraints in classical mechanics, and therefore the normality of the constraint force to the constraint manifold, excludes phenomena as friction.

Combining ( 2) and ( 1) gives the variational equality

(mr -F ) T δr = 0 ∀δr ∈ T K (r), (4)
which is usually referred to as the classical principle of d'Alembert-Lagrange. We now replace the bilateral constraint by unilateral holonomic (geometric) scleronomic constraints g i (r) ≥ 0, i = 1 . . . m, see Figure 2. The position of the point-mass is in this case constrained to the set K = {r ∈ R 3 | g i (r) ≥ 0, i = 1 . . . m}. In the following, we will confine ourselves to the case that K is a convex set. Virtual displacements δr are admissible with respect to the unilateral constraint if δr ∈ T K (r), where T K (r) is now the tangent cone. We define a unilateral constraint as perfect if the constraint force satisfies the normality condition (3)

-R ∈ N K (r), (5)
where N K (r) is the normal cone to the convex set K.

The tangent cone T K (r) and normal cone N K (r) are polar in the sense that for all R and δr satisfying

-R ∈ N K (r), δr ∈ T K (r) (6) it holds that -R T δr ≤ 0. (7)
Hence, for perfect unilateral constraints it holds that the virtual work of the constraint force is nonnegative for admissible virtual displacements, i.e.

R T δr ≥ 0 ∀δr ∈ T K (r). (8)
Combining ( 2) and (1) gives the variational inequality

(mr -F ) T δr ≥ 0 ∀δr ∈ T K (r), (9)
which we will refer to as the principle of d'Alembert-Lagrange in inequality form.

3 The Principle of Hamilton for Non-Impulsive Motion Consider a system S of which we can address each mass-element dm by its position vector ξ ∈ R 3 . The mass-element dm is subjected to external and internal forces dF , which consist of elastic forces, gravitational forces and bilateral holonomic scleronomic constraint forces. Furthermore, the mass-element dm is subjected to unilateral constraint forces dR which impose the unilateral holonomic scleronomic constraint g ≥ 0.

The gap function g depends on the state of the system S and therefore on the position vector ξ ∈ R 3 of all mass-elements of S. The unilateral constraint g ≥ 0 therefore constrains each mass-element addressed by ξ to a set K ξ . The principle of virtual work states that if the virtual work

δW = S δξ T ( ξdm -dF -dR) (10)
vanishes for all virtual variations δξ, i.e.

δW = 0 ∀δξ, (11)
then the system is in dynamic equilibrium. The contact forces dR are assumed to be perfect and it therefore holds that

dR T δξ ≥ 0 ∀δξ ∈ T K ξ (ξ), (12)
where δξ ∈ T K ξ (ξ) denotes a virtual displacement of ξ which is admissible with respect to the unilateral constraint g ≥ 0. The principle of d'Alembert-Lagrange in inequality form (9) for the system S therefore reads as

S δξ T ( ξdm -dF ) ≥ 0 ∀δξ ∈ T K ξ (ξ). ( 13
)
We now choose a set of generalised coordinates q ∈ R n which uniquely describes the state of the system and which is a minimal set of coordinates with respect to all bilateral constraints. Each particle coordinate ξ(q) is therefore completely determined by q. An arbitrary variation δq causes a variation δξ = ∂ξ ∂q δq which is admissible with respect to all bilateral constraints. The unilateral constraint g ≥ 0 constrains the generalised coordinates q to a set K = {q ∈ R n | g ≥ 0} which we assume to be convex. A virtual displacement δq ∈ T K (q) induces a virtual displacement of the particle δξ = ∂ξ ∂q δq ∈ T K ξ (ξ) which is admissible with respect to the unilateral constraint g ≥ 0.

The virtual work of the inertia forces can be rewritten as is the generalised momentum and δT is the variation of the kinetic energy

T = S 1 2 ξT dm ξ = 1 2 qT M (q) q, ( 17
)
with the mass matrix

M (q) = S ∂ξ ∂q T dm ∂ξ ∂q . (18)
Moreover, we introduce

f = S ∂ξ ∂q T dF ( 19
)
as the generalised force. Substitution of ( 15) and ( 19) in the principle of d'Alembert-Lagrange (13) gives

d dt δq T p -δT -δq T f ≥ 0 ∀δq ∈ T K (q), (20)
which is the inequality form of the well known Lagrange central equation [START_REF] Hamel | Elementare Mechanik[END_REF][START_REF] Papastavridis | Analytical Mechanics: A Comprehensive Treatise on the Dynamics of Constrained Systems: for Engineers[END_REF]. The Lagrange central equation ( 20) holds at each non-impulsive time instance t for which the generalised velocities q exist. Hence, we can integrate the central equation ( 20) over a non-impulsive timeinterval I = [t 0 , t f ] which gives

δq T p t f t0 - I δT + δq T f dt ≥ 0 ∀δq ∈ T K (q).
(21) Moreover, we assume the generalised force f to be a potential force -f = ∇V (q) where V (q) is the potential energy and δV = -δq T f . Hence, by defining the Lagrange function L = T -V we arrive at

δq T p t f t0 - I δLdt ≥ 0 ∀δq ∈ T K (q). ( 22
)
If the boundary conditions are fixed, then the variation δq(t) vanishes at t = t 0 and t = t f and we are allowed to interchanging the order of integration and variation such that

-δ I Ldt ≥ 0 ∀δq ∈ T K (q), (23)
with the boundary conditions q(t 0 ) = q 0 , q(t f ) = q f , which is the principle of Hamilton in inequality form for a non-impulsive time-interval I = [t 0 , t f ]. From the principle of Hamilton we can derive the Euler-Lagrange equations in inequality form by evaluating the variation in (23) as

- I ∂L ∂q - d dt ∂ q δqdt ≥ 0 ∀δq ∈ T K (q).
(24) Hence, for each time-instance t ∈ I the variational inequality

- ∂L ∂q - d dt ∂L ∂ q δq ≥ 0 ∀δq ∈ T K (q) (25)
should hold, which can be cast into the differential inclusion

M (q) q -h(q, q) -f R = 0, -f R ∈ N K (q), (26)
where the mass matrix M (q) is defined by ( 17) and the vector

h = - d dt M (q) q + ∂T ∂q - ∂V ∂q T ( 27
)
contains all smooth forces.

A Weak Principle of Hamilton

In this section we derive a weak principle of Hamilton in inequality form for impulsive motion by directly incorporating the impulsive dynamics in the principle of virtual work. Concepts of measure and integration theory appear to be very useful in this respect. We assume that the position ξ(t) is an absolutely continuous function in time and that the velocity ν(t) of the mass element is a function of locally bounded variation without singular terms (which is sometimes called a function of 'special bounded variation'). This implies the following:

1. At each time-instance t we can define a left and right velocity

ν + (t) = lim τ ↓0 ξ(t + τ ) -ξ(t) τ , ν -(t) = lim τ ↑0 ξ(t + τ ) -ξ(t) τ .
(28)

2. The differential measure dξ of the position ξ(t) contains only a density ν with respect to the Lebesgue measure dt dξ = νdt.

(29)

For almost all t we can define a velocity ν(t) = ξ(t).

3. The differential measure dν of the velocity ν(t) contains a density with respect to the Lebesgue measure dt and with respect to the atomic measure dη, i.e.

dν = νdt + (ν + -ν -)dη. ( 30
)
Let as before dR be the unilateral contact force, associated with a unilateral holonomic scleronomic constraint, and dF be the external force on a mass element dm with position ξ. The non-impulsive dynamics of this mass element is therefore described by the equation of motion dm ξ -dF -dR = 0.

(31)

The impulsive dynamics is described by the impact equation

dm(ν + -ν -) -dP = 0, (32)
where dP is the unilateral impulsive force on the mass element dm. The introduction of the differential measure dν allows us to combine the equation of motion (31) and the impact equation ( 32) in a single formula

dm dν -dF dt -dR dt -dP dη = 0, (33)
which is an equality of measures and which should be understood in the sense of integration. The measure dν, which has a density with respect to the atomic measure, is by definition a mapping on the space of continuous functions, i.e. I f T dν only makes sense if f ∈ C 0 (I, R 3 ). The equality of measures immediately leads to a principle of virtual work in differential measures

S δξ T (dm dν -dF dt -dR dt -dP dη) = 0 (34) for all δξ ∈ D 1 (I, R 3 ). The class D 1 (Ω, R) = {y ∈ AC(Ω, R), y ′ ∈ BV (Ω, R
)} comprises all functions on the domain Ω which are absolutely continuous and piecewise C 1 as well as absolutely continuous functions with an accumulation point. We put ξ(ε, t) = ξ(t) + εω(t) with ω ∈ D 1 (I, R 3 ). Often, we will suppress the explicit notation of the time-dependence and simply write ξ(ε) = ξ + εω. For this family of test functions it holds that δξ = ξε (0)δε = ωδε is continuous in time wheras δν = ωδε is not. Moreover, note that ξ(ε) converges for ε → 0 to ξ in the weak norm • 1 , with

y 1 = max x∈Ω |y(x)| + ess sup x∈Ω |y ′ (x)|. (35)
We now make the assumption that both the contact force dR and the contact impulse dP are perfect unilateral constraint forces/impulses, i.e.

-dR

N K ξ (ξ), -dP ∈ N K ξ (ξ). ( 36
)
The virtual work of the contact force dR and the contact impulse dP is therefore non-negative for admissible virtual displacements

δξ T dR ≥ 0, δξ T dP ≥ 0 ∀δξ ∈ T K ξ (ξ). ( 37
)
This bring us to the principle of d'Alembert-Lagrange in inequality form for differential measures

S δξ T (dm dν -dF dt) ≥ 0 ∀δξ = ωδε ∈ T K ξ (ξ) (38
) in which we explicitly write that the virtual displacements δξ are of the form ωδε. A mass element has a constant mass dm. Taking the differential measure-intime of the term δξ T dm ν and applying the chain rule gives d(δξ T dm ν) = δν T dm 1 2 (ν + + ν -)dt + δξ T dm dν, (39) in which we used the equality d(δξ) = dωδε = ωdtδε = δν dt. Moreover, because ν(t) = 1 2 (ν + (t) + ν -(t)) for almost all t, it holds that 1 (ν + + ν -)dt = νdt for Lebesgue integration. We therefore arrive at the variational inequality as being the variation of the kinetic energy T . However, the variations δν are not arbitrary as ν(t) = ν(t) + δν(t) is of the form ν(t) = ν(t) + ε ω(t). The variation δT therefore reduces to the Gâteaux derivative dT (ν; δν) = ∇T δν. Similar as before, we choose generalised coordinates q(t) which form a minimal set of coordinates with respect to the bilateral constraints and with which we can uniquely describe the position ξ(q) of each mass element dm. Moreover, we introduce generalised velocities u(t), which are assumed to be of locally bounded variation, and which are such that dq = udt. The kinetic energy T is a function of ν(q, u) and we can therefore write T as a function T (q, u). Hence, it holds that δT = T q δq + T u δu. The variations δq and δu are not totally arbitrary as they are of the form δq(t) = q(t) + εw(t) and δu(t) = u(t) + ε ẇ(t), where ω = ∂ξ/∂q w. Using the generalised momentum

p = S ∂ξ ∂q T νdm (42)
and the generalised force ( 19) we transform the principle of d'Alembert-Lagrange into

d δq T p -δT dt-δq T f dt ≥ 0 ∀δq = wδε ∈ T K (q), (43)
which is the Lagrange central equation in differential measures. As before, we integrate over a time-interval I = [t 0 , t f ] and consider the generalised force f = -∇V (q) to be a potential force, which yields

δq T p t f t0 - I δLdt ≥ 0 ∀δq = wδε ∈ T K (q), (44)
where L = T -V . Finally, taking fixed boundary conditions at t 0 and t f we obtain

-δ I Ldt ≥ 0 ∀δq = wδε ∈ T K (q). (45)
As mentioned before, the variations δq are not totally arbitrary and the variation δs of the action in (45) reduces to the Gâteaux derivative ds(q; δq). Consequently, we arrive at a weak form of the principle of Hamilton in inequality form

-ds(q; δq) ≥ 0 ∀δq = wδε ∈ T K (q), (46)
with the action integral

s(q) = I Ldt. ( 47
)
This form of the principle of Hamilton is the condition of a weak local extremal of the action s(q) with the weak norm • 1 .

Conclusions

In this paper we derived a weak form of Hamilton's principle as a variational inequality. There are two interesting things to remark at this derivation. First of all, by making use of differential measures we are able to treat the impulsive and non-impulsive dynamics simultaneously. This means that we do not need to use the usual Weierstrass-Erdmann corner conditions for broken extremals (see [START_REF] Troutman | Variational Calculus and Optimal Control: Optimization with Elementary Convexity[END_REF]) to treat an impact as has been done in [START_REF] Panagiotopoulos | Analytical mechanics: Addendum I. Inequality constraints with elastic impacts. The convex case[END_REF][START_REF] Panagiotopoulos | Inequality constraints with elastic impacts in deformable bodies. The convex case[END_REF]. Secondly, as we want to use the principle of virtual work for differential measures (34) we are forced to consider test functions ξ(ε, = ξ(t) + εω(t) of which the variation δξ = ωδε is time-continuous. A test function of this form can have kinks (and therefore impacts), but this class of test functions does not include a family of curves which only varies the impact time, i.e. it is weak. Accordingly, for this weak form of the principle of Hamilton we only have to satisfy the first Weierstrass-Erdmann corner condition. Hence, we are able to prove the validity of ( 46) by only assuming that the unilateral holonomic constraint is perfect which is expressed by The assumption of energy conservation during the impact, i.e. T + = T -, is therefore not a necessary condition for the validity of (46).

Figure 1 .

 1 Figure 1. Point mass in contact with a bilateral constraint.

Figure 2 .

 2 Figure 2. Point mass in contact with a unilateral constraint.

  ∀δξ = ωδε ∈ T K ξ (ξ). We recognise the second termδT = S δν T dm ν, T = S 1 2 ν T dmν (41)


Etta Grover-Silva
email: etta.grover-silva@mines-paristech.fr

Robin Girard

George Kariniotakis

Multi-temporal Optimal Power Flow for Assessing the Renewable Generation Hosting Capacity of an Active Distribution System

The detailed modeling of distribution grids is expected to be critical to understand the current functionality limits and necessary retrofits to satisfy integration of massive amounts of distributed generation, energy storage devices and the electric consumption demand of the future. Due to the highly dimensional non-convex characteristics of the power flow equations, convex relaxations have been used to ensure an efficient calculation time. However, these relaxations have been proven to be inexact during periods of high RES injection. In this paper additional linear constraints were introduced in the power flow formulation to guaranty an exact relaxation. This convex relaxation is then applied within a multi-temporal algorithm in order to evaluate the benefits of storage grid integration. The case study of a French medium voltage feeder is studied to evaluate the maximum capacity of the grid to host RES sources and the advantages of storage systems in reducing curtailment of RES.

I. INTRODUCTION

The electric distribution system is critical for energy security and economic stability. With the exploration of new energy production solutions including renewable energy sources (RES) and associated storage devices, the architecture and functionality of the current distribution grid has been a subject of high interest. For example, a renewable energy source could include wind turbines, hydro turbines or photovoltaic panels. The simulation of current functionality limits with new decentralized renewable energy generation can in turn indicate the advantages of the automation and control of components in the distribution grid. This automation and control is typically labeled as an active distribution system. The added control possibilities of an active distribution grid can allow for a high capacity of decentralized generators to be installed without major infrastructure upgrades. The capacity of a grid to integrate decentralized generators without violating operational limits can be called a distribution grid hosting capacity. In order to quantify the hosting capacity of a distribution grid, detailed AC power flow models are necessary. However, for AC power flow algorithms to be useful and applicable, computation time must be optimized.

There exists many different techniques of AC power flow modeling. However, power flow calculations are non-linear, non-convex and highly dimensional, which can be extremely computationally intensive. Existing power flow algorithms include the forward/backward sweep, Newton Raphson method [START_REF] Vinkovic | Current-based models of {FACTS} devices for three-phase load-flow calculations using the newtonraphson method[END_REF], [START_REF] Penido | A new tool for multiphase electrical systems analysis based on current injection method[END_REF], fast-decoupled load-flow method [START_REF] Abdel-Akher | Improved three-phase powerflow methods using sequence components[END_REF], [START_REF] Yang | Power flow calculation for unbalanced three-phase distribution network with dgs based on phase-sequence hybrid modeling[END_REF], z-bus matrix construction method [START_REF] Farag | A generalized power flow analysis for distribution systems with high penetration of distributed generation[END_REF], and loop impedance method [START_REF] Wu | A three-phase power flow algorithm for distribution system power flow based on loop-analysis method[END_REF], [START_REF] Chen | Loop frame of reference based three-phase power flow for unbalanced radial distribution systems[END_REF]. The quantification of the current hosting capacity of the electric grid is evaluated in [START_REF] Etherden | Increasing the hosting capacity of distribution networks by curtailment of renewable energy resources[END_REF]. As noted within this study, possible strategies to increase this current hosting capacity include curtailment and dynamic line rating. Other strategies well explored in the literature include the use of storage elements [START_REF] Abdelouadoud | Planning-oriented yearly simulation of energy storage operation in distribution system for profit maximization, voltage regulation and reserve provisionning[END_REF].

The development of smart grid solutions and grid automation has increased the passive hosting capacity by controlling certain components during critical times. However, management algorithms for optimal control must resolve this highly dimensional non-convex power flow problem. Convex optimum power flow relaxation algorithms have been developed to optimize the controllable components while ensuring algorithmic efficiency as seen in [START_REF] Gan | Convex relaxations and linear approximation for optimal power flow in multiphase radial networks[END_REF], [START_REF] Madani | Convex relaxation for optimal power flow problem: Mesh networks[END_REF]. Heuristic methods have also been explored to solve the non-convex power flow equations as seen in [START_REF] Ghasemi | An improved teachinglearning-based optimization algorithm using lvy mutation strategy for non-smooth optimal power flow[END_REF]. However, heuristic algorithms often require a larger calculation burden as noted in [START_REF] Paudyal | Three-phase distribution opf in smart grids: Optimality versus computational burden[END_REF] when compared to convex relaxation algorithms. The family of convex relaxation algorithms that is most commonly used for distribution grids is called a Quadratically Constrained Quadratic Program (QCQP). In multiple studies, the nonconvex power flow equations have been cast as a QCQP as shown in [START_REF] Bai | Semi-definite programming-based method for security-constrained unit commitment with operational and optimal power flow constraints[END_REF], [START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF]. Within the QCQP family, two convex relaxation algorithms exist including the Second-Order Cone Program (SOCP) or the Semi-Definite Program (SDP). An SDP convex relaxation has been proved to be exact under certain conditions by [START_REF] Bose | Equivalent relaxations of optimal power flow[END_REF]. While an SOCP relaxation has also proved to be exact under certain conditions as stated in [START_REF] Farivar | Optimal inverter var control in distribution systems with high pv penetration[END_REF], [START_REF] Gan | Convex relaxations and linear approximation for optimal power flow in multiphase radial networks[END_REF], [START_REF] Li | Exact convex relaxation of opf for radial networks using branch flow model[END_REF]. However, these relaxations have been proven to be inexact during periods of high renewable energy production feeding into the grid due to elevated line losses [START_REF] Abdelouadoud | Iterative linear cuts strenghtening the second-order cone relaxation of the distribution system optimal power flow problem[END_REF]. In order to determine the hosting capacity of a distribution grid, it is critical to have a precise and accurate calculation methodology when RE production is high. An example of this difficulty could be high photovoltaic (PV) production during the summer season. In order to overcome the challenge of inaccurate results at periods of high PV injection, [START_REF] Abdelouadoud | Iterative linear cuts strenghtening the second-order cone relaxation of the distribution system optimal power flow problem[END_REF] presents an AC optimum power flow algorithm that integrates linear cuts implemented in a iterative fashion to ensure an exact and feasible relaxation of the power flow equations. This single phase AC optimum power flow algorithm minimizes losses in the system, while also minimizing the distance between the curtailed production and the real PV production capacity. This methodology allows for the instantaneous assessment of grid operational limits with a certain PV injection and a certain consumption. However, this algorithm does not allow for the assessment of time dependent components such as storage devices. In contrast to the literature, the consideration of multiple time steps permits to properly assess the impact of storage since decisions for charging/discharging involve temporal dependencies. Here, the second order cone program (SOCP) convex relaxation algorithm will be implemented using the same linear cuts in a multi-temporal application. The advantage of a multi-temporal problem formulation is the optimization of the charging and discharging schedule of a given time period, here between 12-96 consecutive time-steps.

The importance of multi-temporal coupling in an optimum power flow algorithm is directly related to the evaluation of the role of storage within an active distribution network. An alternative solution to curtailment on the RES generation is the storage of this energy in order to mitigate grid constraint violations. This stored energy can then be later used to minimize the high consumption peak periods. This multi-temporal coupling is critical in order to evaluate the technical constraints of the storage elements and the possible benefits. In this paper, an example active medium voltage distribution grid will be modeled in order to understand the existing RES hosting capacity with curtailment in comparison with centralized and decentralized storage options. This multi-temporal algorithm will allow the quantification of storage element advantages and their optimum placement.

II. PROBLEM FORMULATION A. Power Flow Model

Two different sets of equations can be used within a power flow model, the bus injection model (BIM) or the branch flow model (BFM). While both models can be effective for various applications, the BFM system of equations will be used due to better convergence characteristics as explained in [START_REF] Gan | Convex relaxations and linear approximation for optimal power flow in multiphase radial networks[END_REF] specifically in relation to a radial network topology. Let the labeled nodes i, j and k be oriented as described in the figure below.

node i q branch ij node j q branch jk node k q - -
When considering a radial distribution system with decentralized PV generation and decentralized storage elements, the power flow equations can be written as shown below.

P ij = P load,j + n k=1 P jk + r ij I 2 ij + P pv,j + P st,j (1)
Q ij = Q load,j + n k=1 Q jk + x ij I 2 ij + Q pv,j + Q st,j (2)
Equation ( 1) and ( 2) describe the balance of power from the upstream and downstream branches where P load , P pv , and P st are respectively the instantaneous consumption, PV production and battery storage charge or discharge at a given timestep. P jk , r jk , Q jk ,x jk and I 2 ij are respectively the power, the resistance, the reactive power, the reactance and the current associated with the branch jk. Q load , Q pv , and Q st are respectively the instantaneous reactive consumption, PV reactive production, and battery reactive power values. The voltage at each node can be calculated by equation [START_REF] Abdel-Akher | Improved three-phase powerflow methods using sequence components[END_REF].

|V j | 2 = |V i | 2 -2(r ij P ij + x ij Q ij ) + (r 2 ij + x 2 ij )I 2 ij ( 3
)
where |V j | is the voltage magnitude at node j. The current of each branch is calculated as shown in equation ( 4).

I 2 ij = P 2 ij + Q 2 ij |V i | 2 (4)

B. System Constraints

The system constraints of an electrical distribution network include maximum and minimum voltage limits, maximum and minimum current limits, and operational limits imposed by individual components. The network voltage and current limits can be described as shown below:

V i ≤ |V i | ≤ V i (5)
where V i and V i are respectively the lower and upper limits of the voltage on a line. Individual components such as the PV and battery systems were modeled through their inverter behavior as described in the following sections. 1) PV Inverter Model: The PV inverter behavior was modeled as active and reactive generation source with upper and lower limits, P pv and P pv respectively on active power and s pv , s pv on total apparent power injection for a given timestep.

P pv ≤ |P pv | ≤ P pv s pv ≤ |s pv | ≤ s pv (6)
where s pv is the total apparent power of the power injection as defined below.

s 2 pv = P 2 pv + Q 2 pv s pv = P 2 pv + Q 2 pv (7)
2) Battery Inverter Model: The battery inverter was modeled as either a power injection or power consumption for a given node at a given timestep. The time coupling variable indicating the state of charge (SOC) was calculated based on a charging and discharging efficiency associated with the inverter.

E batt,t1 = E batt,t0 -η ch P ch,t1 -η dch P dch,t1

P st = P ch + P dch s st ≤ |s st | ≤ s st (8)
where η ch , P ch,t1 , η dch , P dch,t1 , s st and s st is the charging efficiency, power absorbed during charging, discharging efficiency, power injected into the grid, upper limit and lower limit of total apparent power exchange with the grid. P st is the power injection or absorbed by the battery as defined below:

s 2 st = P 2 st + Q 2 st s st = P 2 st + Q 2 st (9)
The power flow equations are non-linear and non-convex. Therefore, when solving a highly dimensional power flow problem, convex relaxation have been used to ensure high performance algorithms.

III. OPTIMUM POWER FLOW FORMULATION

The same SOCP convex relaxation used in [START_REF] Abdelouadoud | Iterative linear cuts strenghtening the second-order cone relaxation of the distribution system optimal power flow problem[END_REF] and [START_REF] Farivar | Optimal inverter var control in distribution systems with high pv penetration[END_REF] is implemented in a multi-temporal model in order to analyze the hosting capacity of a distribution grid. This relaxation entails the relaxation of certain equality constraints and the substitution of certain quadratic terms for linear terms. The equality constraints in question, (4), ( 5) and ( 6), are relaxed ultimately relaxing the magnitude of currents within each branch and using a conic formation on the limitation of active power exchange with the grid. Two new variables are introduced to replace quadratic terms, υ i = |V i | 2 and = |I ij | 2 in order to successfully formulate an SOCP problem as explained in [START_REF] Farivar | Optimal inverter var control in distribution systems with high pv penetration[END_REF].

A. Objective Function

In order to optimize the total system functionality, the objective function is composed of two parts. The first part is the minimization of the total losses of the system. The second part minimizes the curtailment of the PV systems in order to maximize renewable energy consumption within the grid. Therefore the objective function is formulated as seen below:

min n i=1 r ij ij + (P pv ideal,i -P pv i ) (10)

B. SOCP Problem Formulation

The complete SOCP formulation is then found below:

min n i=1 r ij ij + (P pv ideal,i -P pv i ) (11)
s.t.

(6), ( 7), ( 8), ( 9)

P ij = P load j + n k=1 P jk + r ij + P pv j + P st j Q ij = Q load j + n k=1 Q jk + x ij + Q pv j + Q st j υ j = υ i -2(r ij P ij + x ij Q ij ) + (r 2 ij + x 2 ij ) ij ij ≥ P 2 ij + Q 2 ij υ i V 2 i ≤ υ ≤ V 2 i ( 12
)

IV. ALGORITHMIC PERFORMANCE

The optimization of storage elements and PV curtailment was analyzed for a time scale of 0.5 to 4 days. An example urban electric grid in France was studied in order to evaluate the performance of this algorithm. The network studied is a 30 kV 137 nodes typical urban grid topology. Network data was acquired through a partnership with ErDF. Twenty PV systems with equivalent power ratings were placed randomly throughout the grid topology with an associated battery system installed at the same node. The time of execution of the algorithm was recorded for varying time coupling scenarios as shown below.

This multi-temporal coupling allows the optimization of PV curtailment and battery storage utilization for up to a four day period with a satisfactory calculation burden. It also ensures that the relaxation is exact and applies linear cuts to timesteps that are not exact in order to guaranty the exactness of the relaxation.

V. CASE STUDY

An example urban electric grid in France was studied to demonstrate possible algorithm utilization. The grid is composed of 137 nodes at 30 kV and serves as a medium voltage distribution grid feeder. A map of the grid topology can be found below. Consumption load profiles were simulated using an aggregation load simulator as described in [START_REF] Barbier | A novel approach for electric load curve holistic modelling and simulation[END_REF] for each low voltage substation and medium voltage consumer. This load simulator takes into account a mix of residential and commercial customers. Residential consumption is simulated with statistically accurate representations of surface area, electric heating and number of individuals that align with the INSEE building inventory database of France. Industrial load profiles are simulated by assuming a typical industrial activity mix in France from a medium voltage substation. The two medium voltage customers were modeled as office building complexes with typical business hours operation. A peak of 5 MW during the winter and 2.9 MW during the summer was simulated for the medium voltage substation transformer with a maximum apparent power limit of 10 MW. Solar radiation data from a site in the south of France was used to calculate expected PV production as a function of system nominal power.

Four different seasonal scenarios were studied in order to understand a typical annual operation. Within each season, three renewable energy solutions were studied: the hosting capacity of a typical grid topology without storage, the hosting capacity with decentralized storage elements installed at the same nodes as the PV systems and the hosting capacity with centralized storage elements close to the high voltage transformer. The hosting capacity of a feeder was defined as the maximum capacity of PV that does not violate grid constraints without using curtailment. Twenty sites were chosen for PV decentralized installations and nominal power of each system was increased until either current or voltage limits were reached.

VI. RESULTS

A three day simulation period was chosen in order to allow for at least one full cycle of charging and discharging of the storage elements. Three day typical profiles were chosen for four different seasons in order to understand the annual performance of the system. Considering only PV installations without storage elements, the maximum current limit of the lines was reached during low loading periods and high peak PV injection in summer. Multiple capacities for 20 decentralized PV systems were tested until curtailment was necessary to not exceed voltage or apparent power limits of the network. The maximum installed capacity where curtailment was unnecessary was achieved with a PV penetration of 10 MW nominal power installed. A a total of 701 kWh of the 156 MWh produced during a three-day simulation was necessary to ensure maximum apparent power limits were not exceeded. During typical daily profiles for fall, winter and spring, no curtailment was necessary for a 10 MW systems. In all simulations, the summer period was the most critical to monitor for grid stability verification, therefore the rest of the results section will focus on summer production and consumption profiles. A comparison of the net consumption of the high voltage transformer and the total power injected into the high voltage grid from the medium voltage substation during summer periods can be found in figure 2.

The curtailment necessary and the losses on the lines were then compared as seen in figure 3. Added storage capacity was then integrated into the grid in order to quantify the additional hosting capacity possible. Initially storage elements were placed at the same nodes as all PV installations to represent a decentralized storage configuration. Therefore 20 systems of 116 kWh energy storage capacity were modeled through their inverter behavior as described in equation 8. For 10 MW of PV capacity installed, no curtailment was needed within all seasons. The improvement in system performance can be seen in figure 4 as a percentage increase or decrease for each parameter. The same amount of storage capacity was then used in a centralized configuration. This centralized storage was placed very close to the high voltage transformer offering the same services. The necessary curtailement was also reduced to zero for the 10 MW system. However, the centralized battery system also resulted in higher overall line losses. The percentage change due to the presence of a centralized battery system is shown in figure 5. The optimization algorithm proposed in this paper successfully applied a convex relaxation algorithm with linear cuts to a multi-temporal application for battery storage analysis. This algorithm was shown to be effective when studying a distribution system for a 3-4 day time span. The multi-temporal algorithm allows for the assessment of battery storage functionality while taking into account the technology limitations such as charging and discharging efficiency, and maximum injection and absorption. The algorithm also calculates an optimal charging and discharging schedule based on the objective function. A French medium voltage distribution feeder was successfully analyzed to determine the hosting capacity, decentralized and centralized battery systems effects on curtailment. This algorithm could be used for comparison studies between different grid stability control strategies such as real time curtailment, centralized and decentralized battery installations. The use of this algorithm for real time management could also be effective if real time predicted PV production profiles and expected load profiles are used as inputs. The optimization of charging and discharging schedules of batteries can effectively increase hosting capacity of a distribution network by reducing the necessity to curtail PV systems.

Fig. 1 .

 1 Fig. 1. Case study of typical medium voltage electric grid topology

Fig. 2 .Fig. 3 .

 23 Fig. 2. Comparison of net consumption of feeder and total PV injection during a three day period using typical summer profiles

Fig. 4 .

 4 Fig. 4. Percentage increase or decrease due to added storage elemnts at each PV production node. Calculations for a 7 MW and 10 MW systems during a three days period using typical summer profiles

Fig. 5 .

 5 Fig. 5. Percentage increase or decrease due to an added centralized storage elemnt when comparing 7 MW and 10 MW systems during a three days period using typical summer profiles

ACKNOWLEDGMENT

The authors would like to thank ERDF for the provision of some of the data used in this study. Mr. Thibaut Barbier is thanked for his support with the simulated load curves. This work is performed in the frame of the European project Sensible (Grant Agreement 645963) funded in part by the European Commission under the Horizon 2020 Framework Program.


Régis Guinvarc'h

Mohammed Serhir

Fabrice Boust

A compact dual polarized 3:1 bandwidth omnidirectional array of spiral antennas

Keywords: Spiral antenna, wideband array, omnidirectional antennas, compact antenna, radar

A circularly polarized wideband omnidirectional array of spiral antennas is presented. It covers the frequency band from 1.2 GHz up to 3.6 GHz while keeping a good S11 and an omnidirectional coverage. This array achieved a dual polarization using left and right hands archimedean spiral antennas. The gain of the array is enhanced using a unique cavity for the entire array and the realized gain is above the maximum theoretical gain of a single spiral antenna.

I. INTRODUCTION

O MNI directional array can be used in various kind of applications, such as direction finding or passive radar. For most of these applications, it is desirable to get as much information as possible from the various azimuthal space directions. So a wide bandwidth and dual polarization capability are some interesting features. Whether it is for airborne, spaceborne or ground based uses, a compact system is a very good asset.

One way to build wideband arrays is to use wideband elements. Considering dual polarized omnidirectional arrays, this has been studied with horn antennas [START_REF] Wang | Tem horn circular array for wide band pattern notch reconfigurable antenna system[END_REF], bowtie etc.. Spiral antennas have been used in different kind of arrays: linear [START_REF] Serhir | A low-profile cavity-backed dual-polarized spiral antenna array[END_REF], planar [START_REF] Saenz | A dualpolarized wideband planar phased array with spiral antennas[END_REF], conformal [START_REF] Chauvet | Approximated method neglecting coupling for conformal array[END_REF] etc., but not in circular compact array, to the best of our knowledge. Its wide impedance bandwidth, along with its excellent polarization properties and wideband radiation pattern, makes it however an interesting candidate.

This paper describes a circular array of spiral antennas providing an omnidirectional coverage, over a wide frequency bandwidth. The dual polarization capability is ensured through the use of an alternate configuration. The array can emit (or collect) signals in any given direction, but not all the direction at the same time. The focus is put on the compactness so that it can be easily used in many applications. It is associated with a good gain thanks to the design of the array and in particular to the backing cavity. Considering passive applications (so positive gain as the criterion), the diameter of the antenna is only one wavelength at the lowest operating frequency. For active applications (so a reflection coefficient less then -10 dB), a 3:1 bandwidth is guaranteed.

Section II describes the array configuration, along with the array element and the cavity design. Section III provides R. Guinvarc'h is with SONDRA, CentraleSupelec, Universit Paris-Saclay, 91192 Gif-sur-Yvette, France ( regis.guinvarch@centralesupelec.fr) M. Serhir, is with GEEPS, UMR 8507 (CentraleSupelec -CNRS UPSud -UPMC) 11 rue Joliot-Curie 91192 Gif Sur Yvette Cedex -France F. Boust is with SONDRA and ONERA -The French Aerospace Lab, Palaiseau, France Fig. 1. The spiral antennas are in alternate configuration. hcav is 30mm, wcav is 122mm and dcav is 280mm. lsusb is 72.4mm and wsubs is 86mm. measured and simulated results comparison. Finally, in Section IV concluding remarks are presented.

II. ANTENNA DESIGN

When designing a compact omnidirectional array of spiral antennas, there are two typical challenges. At low frequencies, the issue is the size of the antenna. Whatever the criteria (gain, polarization etc.), to reach lower frequencies, a spiral antenna has to get bigger. At the other end of the frequency range, the issue is the beamwidth of the spiral antenna to cover 360 • .

A. An Archimedean spiral antenna

The Archimedean spiral antenna is a well known wideband antenna. Its lower operating frequency f min , given in Eq. 1, is inversely proportional to the antenna diameter D. There has been a lot of work to miniaturize it, however, we will use its simplest version, so without any loading, meandering [START_REF] Lee | Distributed lumped loads and lossy transmission line model for wideband spiral antenna miniaturization and characterization[END_REF] etc.

min = c πD (1)

B. Spiral antennas omnidirectional array

The omnidirectional feature requires N spiral antennas placed regularly over a cylinder and facing different directions. The spiral antennas planes are orthogonal to the radii of the circle, as shown on Fig. 1. The resulting structure presents a rotational symmetry without any edge element. Hence, all the elements have the same properties. As for the design, the main parameter is the number of spiral antenna to be used in the array in order to cover the 360 • over the frequency bandwidth. The main difficulty occurs generally at the highest operating frequency. The spiral antenna with its broad beam was particularly interesting for this purpose. The minimum circumference (r) of the array is directly linked to the size of the spiral element D and so to the lowest frequency of operation f min . The antenna array radius is linked to the number of the array spirals and to the spiral antenna diameter. We can easily infer:

r > N D 2π (2)
For instance, a 5-spiral antenna array implies a 72 • beamwidth to cover 360 • . From simulations of single spiral in free space, Eq. ( 2) leads to a maximum value of 3.4 GHz for a single 72mm spiral, with a minimum radius of 5.7 cm. course, this does not take into account the coupling between the array element and with the cavity. This is studied later with simulation of the full structure.

C. Polarization

Spiral antennas of both left and right hand circular polarization have to be used in alternate configuration to get a dual polarization capability. The number of elements N has then to be doubled to account for this alternate configuration. Consequently, the minimum number of antennas to be used to get dual polarization with the omnidirectional feature is 10. Using a spiral antenna with a diameter of 72 mm, this gives the smallest possible array, with a radius of 11.5 cm.

D. Cavity design and gain improvement

The ground plane has been shaped to follow the decagon shape. Two metallic sides have been added to create a cavity. This ensures an unidirectional beam, with lots of space in the middle of the antenna to put components.

the unidirectional feature, the cavity plays also an important role for the gain enhancement. The gain is generally limited by the size of the antenna aperture, this is particularly important at low frequencies. Using the proposed cavity design, we have created an infinite cavity along the array dimension. It allows the field to spread in a larger space. At high frequency, the current does not reach the end of the arms, so there is no difference. At low frequency, instead of being reflected at the end of the arm, the current can couple onto the next spiral through some kind of transmission line with the ground. S 2,1 reaches -12dB at 1.25 GHz and 1.5GHz. This avoids the energy to be brought back to the generator (which would lead to a bad VSWR), while radiating the wrong polarization; on the contrary, the neighboring spiral being of the opposite polarization and fed at the end of its arm by coupling, it radiates the good polarization. This is illustrated in Fig. 2 with a comparison of the reflection coefficient of the array with a common cavity in dashed purple and with separated cavities (one per antenna) in dashed blue. The S 11 cutoff frequency has been lowered by about 350 MHz. In addition, Fig. 3 shows that the common cavity results in 4 dB gain improvement over this bandwidth, with a lower limit to get a positive gain at 1000 MHz instead of 1250 MHz. As the gain improvement is above 3dB, it means a better matching of the antenna due to the design of the cavity.

From Fig. 3 additional oscillations for the common cavity are noticed on the gain, not on the reflection coefficient. This is due to reflections on the ground plane where it is not parallel to the antenna. It is worth noting that the principle of the common cavity is valid for many arrays made of spiral antennas. In [START_REF] Serhir | A low-profile cavity-backed dual-polarized spiral antenna array[END_REF], we introduced it for simple linear arrays with only 4 elements, the cavity being still of finite size. In [START_REF] Saenz | A dualpolarized wideband planar phased array with spiral antennas[END_REF], we extended the idea to a planar array of circular shape, where the cavity can already be seen as 1-D. We use here this principle, but with the array conformed to a decagone, with the antennas outwardly.

III. MEASUREMENTS RESULTS AND ANALYSIS

We have built a prototype with = 10 elements. This is presented in Fig. 4 and7. The cavity is made from copper and describes the decagone over which a foam material is placed (white part). Its dimensions are taken from [START_REF] Serhir | A low-profile cavity-backed dual-polarized spiral antenna array[END_REF] and are recalled in Fig. 1. This foam material present a permittivity near to 1 and is supposed to not disturb the radiating behavior of the spiral antennas. These 5-turns spiral antennas are printed over FR4 dielectric substrate with r = 4.2 of 1.6 mm thickness. Then the antennas are placed over the foam material. The top view of the antenna array prototype is presented in Fig. 4. In the middle of the photography we see the feeding cables associated to each antenna. Indeed, we used a pair of coaxial cable to excite symmetrically each spiral antenna. Two RG-10 coaxial cables are soldered together and the inner conductor of each cable is connected to the spiral arms. So the impedance between the two inner conductors (so seen from the antenna) is 100 ohms. The input impedance of the spirals when placed over a FR4 substrate with 1.6 mm thickness is 120 Ω as presented in [START_REF] Serhir | A low-profile cavity-backed dual-polarized spiral antenna array[END_REF]. The two opposed-phase signals associated to each coaxial cable are generated using a 180 • Fig. 4. Antenna in the anechoic chamber with an inset of the top view of the antenna, with the feeding cables paired to ensure the matching.

A. Reflection coefficient

The measured and the simulated reflection coefficients are plotted in Fig. 5. It must be mentioned that, thanks to the feeding structure depicted above, the reference impedance is 50 Ω. When not in an array, the reflection coefficient of the spiral antenna should be less than -10 dB starting at 1.32 GHz ( c Dπ ). When placed in the array, the same spiral antenna presents a reflection coefficient S11 below -10 dB after 1.27 GHz. It is satisfactory up to at least 4 GHz. For the passive applications, it is interesting to consider that it may be usable starting around 1 GHz. There is a very good agreement between the measurement and the simulation, especially at low frequencies. The ground plane being at 3 cm from the antenna, the destructive interference should occur around 5 GHz. The active reflection coefficient [START_REF] Holter | On the size requirement for finite phasedarray models[END_REF] is also shown. It is very almost identical to the S11, except at low frequency where the maximum difference is only 2.5 dB.

B. Gain

In Fig. 6 we compare the measured gain of the array and the simulated one including the losses of the feeding system. These losses are quite high (above 4dB) because of the large losses in the hybrid couplers. However, these particular hybrids are interesting as they are quite small and as such, quite easy to include in the structure. Despite these high losses, the system gain is positive, except at the low frequency end. For passive radar at these bands, it is not an issue. There is a deep caused by a beam splitting around 1.5 GHz in simulation orthogonally to the array axis. It is shifted in the measurements to 1.6 GHz and almost canceled (less than 3 dB, see Fig. 9). Along the array axis, the deep does not appear (see Fig. 8). In a first design, we used larger spirals of 78 mm [START_REF] Guinvarc'h | A wideband omnidirectional conformal array for passive radar[END_REF]. The deep was much stronger and not canceled in the measurements. The maximum expected gain 4πS λ 2 is also plotted on Fig. 6, for a surface S of the size of a single antenna substrate in red and for the projected surface of the whole antenna in green. The gain of the antenna is higher than the theoretical maximum gain of a single antenna on half of the bandwidth. It demonstrates the interest of the unique cavity to spread the field on the whole antenna and thus to improve the gain.

C. Radiation patterns

Fig. 8 shows the measured normalized radiation patterns for eight frequencies, from 1200 MHz up to 3900 MHz, and are compared with the simulation results. Simulations and measurements show similar behavior. At broadside, the agreement is very good especially for low frequencies (see Fig. 6). However, out of the main beam there can be some discrepancies. This is probably due to fabrication inaccuracy and measurement uncertainties in the anechoic chamber (seen for higher frequencies), as well as some effect of the supporting structure. This was not taken into account in the simulation as the antenna array is too large to be simulated. In order to fulfill the 360 • coverage requirements, every pattern must have a 360 • /5= 72 • 3-dB beamwidth in the φ = 0 • plane, whatever the frequency. This was achieved for the simulation. But for the prototype, the criteria has to be lowered to a 4-dB beamwidth, which is still a good value given the compact size of the antenna. Fig. 9 plots the measured radiation patterns for four frequencies, orthogonally to the array axis. The beam is quite broad, typical of spiral antenna despite being narrower in this plane cut because of the sides of the cavity. It goes narrower as the frequency increases. At 1.6 GHz, the beam is slightly split. This corresponds to the loss of broadside gain at 1.6 GHz of Fig. 6. So that loss is not due to a mismatching (as could also be seen on Fig. 5). There were no specific requirements over the radiation pattern in this plane.

D. Cross polarization rejection

The cross polarization rejection of Fig. 10 demonstrates a good agreement between simulation and measurement results. The cross polarization rejection is around 10 dB, which is not a satisfactory value for many applications. It is particularly interesting to do the comparison with a single spiral antenna above a ground plane (dashed blue line in Fig. 10). The values are roughly the same, in particular the low frequency behavior is absolutely not affected by the array, despite the conformal shape that could have been expected to alter the polarization. This leads to the conclusion that the array geometry, including the cavity, does not alter the polarization. As a consequence, the polarization could be improved by simply using optimized spirals with slow-wave techniques (see [START_REF] Lee | Distributed lumped loads and lossy transmission line model for wideband spiral antenna miniaturization and characterization[END_REF] for instance).

IV. CONCLUSION

A compact wideband omnidirectional array of spiral antennas has been presented. It is almost one lambda of diameter at the lowest operating frequency. It works from 1.2 GHz up to 3.6 GHz, with a good reflection coefficient and an omnidirectional coverage. The use of a unique cavity has increased the gain above the maximum theoretical gain of one antenna, in addition to put the cavity effects outside the bandwidth. It results in a positive gain from 1 (so below the theoretical cutoff frequency of the spiral) to at least 4 GHz.

Fig. 2 .

 2 Fig. 2. Simulated reflection coefficients for common and separated cavities.

Fig. 3 .

 3 Fig. 3. Simulated gain for common and separated cavities.

Fig. 5 .

 5 Fig. 5. Comparison between measured and simulated reflection coefficients of the prototype (without the matching structure losses). The active reflection coefficient (ARC) is also included.

Fig. 6 .

 6 Fig. 6. Comparison between measured and simulated gain of the prototype (without the matching structure losses).

Fig. 7 .

 7 Fig. 7. Surfaces taken into account to calculate 4πS λ 2 .

Fig. 8 .

 8 Fig. 8. Measured (solid lines) and simulated (dashed lines) radiation patterns for 8 frequencies, along the array axis (xz plane of Fig. 1).

Fig. 9 .

 9 Fig. 9. Measured radiation patterns for 4 frequencies, orthogonal to the array axis (yz plane of Fig. 1).

Fig. 10 .

 10 Fig. 10. Comparison between measured and simulated rejection of cross polarization rejection. The corresponding value for the spiral antenna alone above a ground plan is also indicated for reference. The horizontal dashed line at 15dB indicates an Axial Ratio of 3dB.


G Chen

Bernard Roux

whether they are published or not. The documents may come

Analytical solution and numerical simulation of thermocapillary convection in floating zones

Analytical solution and numerical simulation of thermocapillary convection in floating zones


Heat transfer, thermocapillary convection and interface shapes in floating zone crystal growth with a magnetic field

G. Chen, B.