Mihai Anitescu
email: anitescu@mcs.anl.gov

Alessandro Tasora
email: tasora@ied.unipr.it

An iterative approach for cone complementarity problems for nonsmooth dynamics

Keywords: Iterative methods, cone complementarity problems, LCP, complementarity, contacts, multibody

Aiming at a fast and robust simulation of large multibody systems with contacts and friction, this work presents a novel method for solving large cone complementarity problems by means of a fixed-point iteration. The method is an extension of the Gauss-Seidel and Gauss-Jacobi method with overrelaxation for symmetric convex linear complementarity problems. The method is proved to be convergent under fairly standard assumptions and is shown by our tests to scale well up to 500,000 contact points and more than two millions of unknowns.

Introduction

Mechanisms involving contacts and impacts between parts can be modeled in terms of multibody systems with unilateral constraints. The simulation of rigid contacts entails the solution of nonsmooth equations of motion: the dynamics is nonsmooth because of the discontinuous nature of noninterpenetration, collision, and adhesion constraints [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF].

Devices composed of rigid bodies interacting through frictional contacts are extensively used in many engineering solutions, either featuring a small number of unilateral contacts such as cam-followers and Geneva wheels or including thousands of contacts between a large number of parts, such as in the cases of palleting machines, vibratory feeders, size-segregation devices, CVT chains, and pebble reactors. Robust and efficient simulation software is mandatory, since the proper operation of these devices relies on the formation and loss of frictional contacts that cannot be easily studied by analytical methods. Given the presence of discontinuities, however, a straightforward application of numerical methods for ordinary differential equations is impracticable.

One of the most popular approaches to nonsmooth dynamics is the integratedetect-restart method, which adopts traditional DAE or ODE integration on piecewise integrals (Caratheodory integrals) [START_REF] Haug | Computer Aided Kinematics and Dynamics of Mechanical Systems[END_REF][START_REF] Haug | Dynamic mechanical systems with coulomb friction, stiction, impact and constraint addition-deletion[END_REF]. Nonetheless, this scheme, though reliable for systems with one degree of freedom, may fail when handling multiple unilateral constraints, because there is no way to guarantee an upper bound on the number of subproblems to solve in finite time intervals [START_REF] Trinkle | On dynamic multi-rigid-body contact problems with coulomb friction[END_REF].

Another popular approach is represented by regularization strategies, which model contacts by means of many compliant spring-dashpot linkages. This approach requires little effort on the programming side and allows the adoption of normal ODE or DAE integrators [START_REF] Donald | On the motion of compliantly connected rigid bodies in contact: a system for analyzing designs for assembly[END_REF][START_REF] Song | Analysis of rigid-body dynamic models for simulation of systems with frictional contacts[END_REF][START_REF] Song | A semi-implicit time-stepping model for frictional compliant contact problems[END_REF][START_REF] Pang | Convergence of time-stepping method for initial and boundary-value frictional compliant contact problems[END_REF]. Because of the high stiffness affecting the explicit integration, however, this method may require prohibitively small time steps to achieve alpha stability; moreover, the need to tune additional parameters on case-by-case basis is not welcome by end users.

These issues motivate the investigation of innovative numerical methods that can deal with multiple frictional contacts, even in case of thousands, if not millions, of moving parts. To that end, much attention was drawn by time-stepping approaches that produce weak solutions of the differential variational inequality (DVI) that describes the continuous time motion of rigid bodies with collision, contact, and friction. The DVI as a problem formulation was recently introduced in full generality and classified by differential index [START_REF] Pang | Differential variational inequalities[END_REF][START_REF] Pang | Solution dependence on initial conditions in differential variational inequalities[END_REF], though earlier numerical approaches based on DVI formulations do exist [START_REF] Moreau | Standard inelastic shocks and the dynamics of unilateral constraints[END_REF][START_REF] Marques | Differential Inclusions in Nonsmooth Mechanical Problems: Shocks and Dry Friction, volume 9 of Progress in Nonlinear Differential Equations and Their Applications[END_REF][START_REF] Lotstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]. Recent work on time-stepping schemes has included both acceleration-force linear complementarity problem (LCP) approaches [START_REF] Baraff | Issues in computing contact forces for non-penetrating rigid bodies[END_REF][START_REF] Pang | Complementarity formulations and existence of solutions of dynamic multi-rigid-body contact problems with coulomb friction[END_REF][START_REF] Trinkle | On dynamic multi-rigid-body contact problems with coulomb friction[END_REF] and velocity-impulse LCP-based time-stepping methods [START_REF] Stewart | An implicit time-stepping scheme for rigidbody dynamics with inelastic collisions and Coulomb friction[END_REF][START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF][START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF][START_REF] Anitescu | Time-stepping for three-dimensional rigid-body dynamics[END_REF].

The introduction of inequalities in time-stepping schemes for DVI, coupled with a polyhedral approximation of the friction cone, leads to linear complementarity problems (LCP) [START_REF] Stewart | An implicit time-stepping scheme for rigidbody dynamics with inelastic collisions and Coulomb friction[END_REF], which are systems of complementary inequalities to be satisfied simultaneously [START_REF] Cottle | The Linear Complementarity Problem[END_REF]. These complex LCP problems must be solved at each time step in order to advance the integrator [START_REF] Lotstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF][START_REF] Stewart | An implicit time-stepping scheme for rigidbody dynamics with inelastic collisions and Coulomb friction[END_REF].

If the simulation entails a large number of contacts and rigid bodies as is the case of part feeders, packaging machines, and conveyor belts, the computational burden of classical LCP solvers can be significant. Indeed, a wellknown class of approaches to LCP problems is based on simplex methods, also known as direct or pivoting methods, originating from the algorithms of Lemke and Dantzig [START_REF] Cottle | Complementary pivot theory of mathematical programming[END_REF]. However, these methods may exhibit an exponential worst-case complexity [START_REF] Baraff | Fast contact force computation for nonpenetrating rigid bodies[END_REF]. Our experience shows that, in spite of deep optimizations [START_REF] Tasora | Un nuovo metodo del simplesso per il problema di complementarit lineare mista in sistemi multibody con vincoli unilateri[END_REF], simplex methods still cannot practically handle multibody systems with more than one hundred colliding bodies.

In the three-dimensional case, the Coulomb friction at contact points without the use of a polyhedral approximation leads to a more complex nonlinear complementarity problem (NCP). The use of a polyhedral approximation make possible the use of typical LCP solvers [START_REF] Stewart | An implicit time-stepping scheme for rigidbody dynamics with inelastic collisions and Coulomb friction[END_REF][START_REF] Trinkle | On dynamic multi-rigid-body contact problems with coulomb friction[END_REF][START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF]. Artificial anisotropy, however, affects friction because friction cones become faceted friction pyramids. In addition, such finite approximation of cones results in a far larger problem (insofar as the number of constraints) and has a negative impact on the performance of LCP solvers, which is already critical in general.

To circumvent the difficulties posed by increasing complexity of classical LCP solvers and the increased size and inaccuracy introduced by polyhedral approximation, we have developed a novel solution method, based on a fixed-point iteration with projection on a convex set, that can directly solve large cone complementarity problems with low computational overhead. The method is based on a time-stepping formulation that solves at every step a cone constrained optimization problem [START_REF] Anitescu | Optimization-based simulation of nonsmooth rigid multibody dynamics[END_REF]. The time-stepping scheme was proven to converge in a measure differential inclusion sense, to the solution of the original continuous-time DVI. At every step, we solve the cone complementarity problem (CCP) that results from the optimality conditions of the cone constrained optimization problem. Note that the same formulation has been recently used as the basis of a quasistatic frictional contact model with local compliance [START_REF] Pang | On a continuous-time quasistatic frictional contact model with local compliance[END_REF] that also results in a CCP.

In systems with bilateral constraints only, our method reduces to a stationary Gauss-Seidel or Gauss-Jacobi method with successive over-relaxation [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF]. For the original NCP formulation, Gauss-Seidel methods have been successfully used for thousands of rigid bodies in contact [START_REF] Moreau | Numerical treatment of contact and friction: The contact dynamics method[END_REF][START_REF] Jourdan | A Gauss Seidel like algorithm to solve frictional contract problems[END_REF]. For these methods, however, no convergence theory exists, except for small friction coefficients, whereas our methods converges under certain conditions that do not include a small friction assumption. If the CCP is solved without overrelaxation, then our method shares certain features with the block coordinate descent method with convex constraints on the variables in a block [START_REF] Tseng | A coordinate gradient descent method for nonsmooth separable minimization[END_REF].

Among the most promising applications of this method are dynamical analysis of large scenarios comprising thousands of colliding bodies, as in the case of the simulation of pebble bed nuclear reactors, granular flows, masonry stability analysis, robotics, and CAD/CAM/CAE simulations of complex devices (Fig. 1), which to date are strongly limited by computational complexity issues, even on supercomputers.

Optimization-based Time-Stepping Scheme

In the following, we present our contact model, and we compare it to previous approaches. The object of study is a system of rigid bodies, described by state variables and contact and frictional constraints.

System representation

At a time t, the position of the system is described by generalized coordinates q(t) ∈ R m (which may include rotational coordinates that cannot be defined over a subspace homeomorphic to R n , for some n), and generalized velocities v(t) ∈ R m . In classical mechanics, v(t) is continuous, and we can write dq/dt = Γ (q)v, where Γ (q) is a matrix that connects the generalized velocities to the derivatives of the generalized coordinates.

In three dimensions, the position of a rigid body is described by the position x, y, z of the center of mass and a 3 × 3 orthogonal matrix A ∈ SO(3, R) that represents the rotation of a frame attached to the body with respect to a fixed-world frame. The A matrix could be function of three parameters ∈ R 3 (Cardano angles, Euler angles, etc.), but this may cause numerical problems such as singularities in transformations. To overcome such problems, we adopt four-dimensional unitary quaternions η ∈ S 3 ⊂ H, though their space is not homeomorphic to R 3 .

Here we assume that A can be represented smoothly by three parameters, φ, θ, ζ: this parameterization is valid only locally, while for the global cumulative rotations we use quaternions. This reparameterization does not affect the dynamics [START_REF] Haug | Computer Aided Kinematics and Dynamics of Mechanical Systems[END_REF]. Therefore a system with n bodies in three dimensions is represented by m = 6n coordinates.

Nonpenetration constraints

Two rigid bodies should not penetrate, and, if they are in contact, there should be friction acting at the interface. To enforce the nonpenetration constraint, we assume that there exists a function Φ(q), which we call the gap function, that satisfies

Φ(q) = ⎧ ⎨ ⎩ > 0 if
the bodies are separated, = 0 if the bodies touch each other, < 0 if the bodies are interpenetrating.

(1)

For such a function, the nonpenetration constraint becomes Φ(q) ≥ 0.

An example of such a mapping is the signed distance function [START_REF] Kim | Deep: Dual-space expansion for estimating penetration depth between convex polytopes[END_REF], which is differentiable when the bodies are smooth and convex, at least up to some value of the interpenetration [START_REF] Anitescu | Formulating 3d contact dynamics problems[END_REF]. For most cases, even simple ones involving the relative position of two spheres, a differentiable signed distance function cannot be defined for all values of q. The fact that Φ(q) can be differentiably defined only on a neighborhood of the set Φ(q) ≥ 0 can be accommodated at the cost of making the analysis substantially more involved [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF]. To simplify our discussion, we make the following assumption.

Differentiability of geometrical constraint data assumption:. Any contact is described by a gap function Φ(q) that is everywhere twice continuously differentiable.

Frictional constraints

In this work we describe the frictional constraints by conic constraints, which are an extension of complementarity models discussed in [START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF][START_REF] Stewart | An implicit time-stepping scheme for rigidbody dynamics with inelastic collisions and Coulomb friction[END_REF].

The Coulomb friction model

The model we represent and approximate is the Coulomb friction model. If a position q is feasible and the contact is active, that is, Φ(q) = 0, then at the contact we have a normal force and a tangential force.

Let n be the normal at the contact, pointing toward the exterior of the body, and let t 1 and t 2 be the tangents at the contact. Here n, t 1 , t 2 are mutually orthogonal vectors of length one in three dimensions. The vectors n, t 1 , and t 2 are a function of the position q. In the following we use v, v to refer to velocities, and the subscripts u , v to refer to quantities related to the the two linearly independent tangential directions at a given contact.

The reaction force is impressed on the system by means of multipliers γ n ≥ 0, γ u , and γ v . The normal component of the force is F N = γ n n, and the tangential component of the force is

F T = γ u t 1 + γ v t 2 .
The Coulomb model consists of the following constraints:

γ n ≥ 0, Φ(q) ≥ 0, Φ(q) γ n = 0, μ γ n ≥ γ 2 u + γ 2 v , ||v T || μ γ n -γ 2 u + γ 2 v = 0, F T , v T = -||F T || ||v T || (2)
where v T is the relative tangential velocity at contact. The effect of the friction over the dynamical system is defined by the friction coefficient μ ∈ R + , that typically has a value between 0 and 1 for most materials. 1The first part of the constraint can be restated as

F = F N + F T = γ n n + γ u t 1 + γ v t 2 ∈ K,
where K is a cone in three dimensions, whose slope is arctan(μ).

The constraint F T , v T = -||F T || ||v T || requires that the tangential force be opposite to the tangential velocity. This results in the reaction force being dissipative. In fact, an equivalent convenient way of expressing this constraint is by using the maximum dissipation principle [START_REF] Stewart | An implicit time-stepping scheme for rigidbody dynamics with inelastic collisions and Coulomb friction[END_REF][START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of Painleve's problems[END_REF][START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF] (

γ u , γ v ) = argmin √ b γ 2 u +b γ 2 v ≤μb γ n ( γ u t 1 + γ v t 2 ) T v T .
These constraints are represented by mapping the vectors n, t 1 , t 2 from contact coordinates to generalized coordinates [START_REF] Anitescu | Formulating 3d contact dynamics problems[END_REF].

For example, if we have a two-body system, then the generalized coordinates in the three-dimensional space are embedded in a twelve-dimensional space by using the coordinates

x 1 , y 1 , z 1 , φ 1 , θ 1 , ζ 1 , x 2 , y 2 , z 2 , φ 2 , θ 2 , ζ 2 .
For a three-dimensional vector v, the mapping to generalized coordinates is

v → ⎛ ⎜ ⎝ v r 1 × v -v -r 2 × v ⎞ ⎟ ⎠ ,
where r 1 and r 2 are the relative positions of the contact point with respect to the centers of mass of the two bodies [START_REF] Anitescu | Formulating 3d contact dynamics problems[END_REF]. Using this mapping, we denote the generalized vector version of n, t 1 , t 2 by D n , D u , D v . One unfortunate side effect of generalized coordinates mapping is that, in the new coordinates, D n , D u , D v cease to be mutually orthogonal.

If v is the generalized velocity, the tangential velocity satisfies the following

t T 1 v T = v T D u , t T 2 v T = v T D v .
In generalized coordinates, the Coulomb model thus becomes

F N = γ n D n , F T = γ u D u + γ v D v , γ n ≥ 0, Φ(q) ≥ 0, γ n Φ(q) = 0, μ γ n ≥ γ 2 u + γ 2 v , argmin √ b γ 2 u +b γ 2 v ≤μb γ n ( γ u D u + γ v D v ) T v = ( γ u , γ v ) .

The overall dynamical model

The other dynamical data needed for the model are the mass matrix M (q), the external force f e (t, q, v), and the term f c (q, v) which contains the effect of centrifugal and Coriolis forces. The mapping f c (q, v) is continuously differentiable and satisfies [START_REF] Anitescu | Time-stepping schemes for stiff multi-rigidbody dynamics with contact and friction[END_REF] v T f c (q, v) = 0 ∀q, v.

This equation implies that Coriolis and centrifugal forces do not provide any net work to the rigid multibody system. To simplify notation, we also make the following assumption.

Constant mass matrix assumption: The mass matrix M (q) ∈ R m×m is positive definite and constant. This assumption is satisfied in two dimensions and three dimensions if we use the Newton-Euler formulation in body coordinates [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF].

With this definition, we can define the total force

f t (t, q, v) = f e (t, q, v) + f c (q, v). (3)
Assume now that we have p potential contact constraints, which are enforced by the nonpenetration constraints Φ i (q) ≥ 0, i = 1, 2, . . . , p.

In the following, we denote by the superscript i the data associated to the potential contact i. The continuous model is the following differential variational inequality [START_REF] Pang | Differential variational inequalities[END_REF]:

M dv dt = i=1,2,...,p γ i n D i n + γ i u D i u + γ i v D i v + f t (t, q, v) dq dt = Γ (q) v γ i n ≥ 0 ⊥ Φ i (q) ≥ 0, i = 1, 2, . . . , p γ i u , γ i v = argmin μ i b γ i n ≥ √ (b γ i u ) 2 +(b γ i v ) 2 γ u D i u + γ v D i v T v, i = 1, 2, . . . , p.
(4) Here Γ (q) is used to transform the generalized velocities into derivatives of the generalized positions. For instance, when dealing with rotations, Γ (q) can be a linear mapping from three-dimensional angular speeds into fourdimensional time derivatives of unit quaternions.

The Coulomb model used in this work is the predominant model used in the engineering literature to describe dry friction. Unfortunately, the model may be inconsistent: there exist configurations for which the model does not have a solution [START_REF] Baraff | Issues in computing contact forces for non-penetrating rigid bodies[END_REF][START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF]. This situation has led to the need to explore weaker formulations where the forces are measures and Newton's law is satisfied in a measure differential inclusion sense [START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF]. It has been shown that solutions in that sense do exist and can be found by time stepping schemes [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of Painleve's problems[END_REF].

We will consider all collisions that appear during the simulation of the inelastic type. Therefore, they are naturally treated by the time-stepping scheme through a change of active set without the need to modify the algebraic expression of the scheme.

Time stepping scheme

We now define a stepping scheme for the continuous time formulation. We start at the time t (l) , position q (l) , and velocity v (l) with time step h. The scheme is expressed by the following equation problem with equilibrium constraints:

M (v (l+1) -v l ) = j∈A(q (l) , ) γ i n D i n + γ i u D i u + γ i v D i v +hf t (t (l) , q (l) , v (l) ) ( 5 ) 0 ≤ 1 h Φ i (q (l) ) + ∇Φ i T v (l+1) (6)
⊥ γ i n ≥ 0, i ∈ A(q (l) , ) γ i u , γ i v = argmin μ i γ i n ≥ √ (γ i u ) 2 +(γ i v ) 2 v T (γ u D i u + γ v D i v ) (7)
i ∈ A(q (l) , )

q (l+1) = q (l) + hΓ (q (l) )v (l+1) , (8)
where

A(q, ) = i i ∈ {1, 2, . . . , p} , Φ i (q) ≤ . ( 9
)
We have denoted by γ s the constraint impulses of a contact constraint, that is,

γ s = h γ s , for s = n, u, v. The 1 h Φ i (q (l)
) term achieves constraint stabilization and its effect is amply discussed in [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF].

In previous work, we have shown that the scheme is convergent, as the time step h goes to 0 to the solution of a measure differential inclusion [START_REF] Anitescu | Optimization-based simulation of nonsmooth rigid multibody dynamics[END_REF]. Solutions of the subproblems, when the nonlinear constraint is approximated by a piecewise linear cone, can be found by Lemke's algorithm [START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF]. Nonetheless, in [START_REF] Anitescu | A fixed-point iteration approach for multibody dynamics with contact and friction[END_REF] we have also demonstrated that, as the number of constraints in the problem increases, the computational cost of Lemke's method increases far faster than linearly with the size of the problem. As an alternative we proposed to solve the problem as a monotone optimization by introducing a relaxation over the complementarity constraints, that is, we modified the time-stepping scheme by replacing the equation ( 6) with:

0 ≤ 1 h Φ i (q (l) ) + ∇Φ i T v (l+1) -μ i (D i,T u v) 2 + (D i,T v v) 2 ⊥ γ i n ≥ 0, i ∈ A(q (l) , ). (10)
We note that the modified formulation does fit the paradigm put forth in [START_REF] Pang | A unified approach to frictional contact problems[END_REF].

We also note, however, that, at least in its most evident obvious formulations the problem would violate either assumption (H2) or assumption (H4) that are used in that reference to prove existence of solutions, though it is also clear that even mild regularizations would satisfy those assumptions.

In this work, we do not discuss elastic or partially inelastic collision. This is equivalent to considering a 0 restitution coefficient. It is conceivable that the approach in [START_REF] Potra | A linearly implicit trapezoidal method for integrating stiff multibody dynamics with contact and friction[END_REF][START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF] can be adapted to our scheme if restitution is needed, but that is an issue that needs further study, due to the effective normal compliance that this scheme presents, as discussed below. What ever the modification, however, it is likely to result in a cone complementarity onestep problem whose matrix has similar structure to the one considered here. How to solve that problem, as opposed to the time stepping scheme itself, is the focus of this article.

We have shown in [START_REF] Anitescu | Optimization-based simulation of nonsmooth rigid multibody dynamics[END_REF] that, as h → 0, the solution of the modified timestepping scheme will approach the solution of the same measure differential inclusion as the original scheme. In addition, we have shown that the iterates produced by the modified scheme approach the ones of the original scheme

provided that μ i γ i n (D i,T u v) 2 + (D i,T v v) 2 << 1 [4]
. We remark that this regime is precisely the one in which pebble bed simulators operate [START_REF] Rycroft | Analysis of granular flow in a pebble-bed nuclear reactor[END_REF], an application that motivates our second example, as well as other granular flow applications.

The physical meaning of the modified scheme was described to some extent in [START_REF] Anitescu | Optimization-based simulation of nonsmooth rigid multibody dynamics[END_REF], and is connected to the microscopic description of contact with friction as rigid piecewise linear asperities whose tangent of the side angle is equal to the macroscopic friction coefficient. The effect of the modification on the iterates of the time stepping scheme is that it allows for some amount of normal motion at the contact even when the body is supposed to be in contact. In some sense, the scheme allows for a "boundary layer", effectively normal compliance, at the contact whose size is proportional to μv T h (but that otherwise does not depend on any parameter, as penalty schemes do). This interpretation can be demonstrated on a two-dimensional example assuming persistent contact that is already established. Of course, the modification [START_REF] Cottle | Complementary pivot theory of mathematical programming[END_REF] may not be suitable for all simulations with contact and friction and a good example of an unsuitable configuration is the case where two bodies are in initial contact with large tangential velocity, as is discussed in [START_REF] Anitescu | Optimization-based simulation of nonsmooth rigid multibody dynamics[END_REF].

For the rest of the paper we use the modified scheme, which uses (10) instead of (6).

Cone complementarity formulation

If we now write the optimality conditions for the equilibrium constraint in [START_REF] Anitescu | Time-stepping for three-dimensional rigid-body dynamics[END_REF], we obtain that, for any i ∈ A(q (l) , ), there exists a Lagrange multiplier λ i such that

λ i γ i u = -D i,T u v, λ i γ i v = -D i,T v v, λ i ≥ 0 ⊥ μ i γ i n -(γ i u ) 2 + (γ i v ) 2 ≥ 0. ( 11
)
The first two equations imply that

λ i (γ i u ) 2 + (γ i v ) 2 = D i,T u v 2 + D i,T v v 2 ,
while the last equation implies that

0 = λ i (γ i u ) 2 + (γ i v ) 2 μ i γ i n -(γ i u ) 2 + (γ i v ) 2
and, in turn, that

μ i γ i n D i,T u v 2 + D i,T v v 2 = λ i γ i u 2 + γ i v 2 . ( 12
)
We now define, for i ∈ A(q l , ) the vectors

u i = (u i 1 , u i 2 , u i 3 ), w i = (γ i n , γ i u , γ i v ) u i 1 = 1 h Φ i (q (l) ) + ∇Φ i T v (l+1) , u i 2 = D i,T u v, u i 3 = D i,T v v.
We calculate the scalar product

u i T w i = γ i n 1 h Φ i (q (l) ) + ∇Φ i T v (l+1) + γ i u D i,T u v + γ i v D i,T v v (10),(11) = μ i γ i n D i,T u v 2 + D i,T v v 2 -λ i γ i u 2 + γ i v 2 = 0,
which implies that u i T w i = 0, and thus u i ⊥ w i . (13) We now define the cones

Λ i = (x, y, z) ∈ R 3 |x ≥ μ i y 2 + z 2 , FC i = (x, y, z) ∈ R 3 |μ i x ≥ y 2 + z 2
It immediately follows that Λ i is the negative polar cone of FC i , that is, ũ ∈ Λ i and w ∈ FC i imply that ũT w ≥ 0. Then ( 7), [START_REF] Cottle | Complementary pivot theory of mathematical programming[END_REF], and (13) imply that the following set of cone complementarity constraint holds:

-u i ∈ FC i• ⊥ w i ∈ FC i , i ∈ A(q (l) , ), (14)
where we denote by C • the polar cone of a given cone C, that is,

C • = { x ∈ R m | x, y ≤ 0, ∀y ∈ C}.
We now define the vector k(l) = M v (l) + hf t (t (l) , q (l) , v (l) ).

Then, equations ( 15) and ( 14), together with [START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF] and the definition of the vectors u i and w i , result in the following problem:

M v (l+1) = k(l) + i∈A(q (l) , ) γ i n D i n + γ i u D i u + γ i v D i v , i ∈ A(q (l) , ) 1 h Φ i (q (l) ) + ∇Φ i T v (l+1) , D i,T u v (l+1) , D i,T v v (l+1) ∈ -FC i• ⊥ (γ i n , γ i u , γ i v ) ∈ FC i . (16)
We denote by n A the number of elements in the set A(q l , ). We then define the following vectors:

b ∈ R 3n A = 1 h Φ i 1 (q (l) ), 0, 0, 1 h Φ i 2 (q (l) ), 0, 0, . . . , 1 h Φ i n A (q (l) ), 0, 0 d ∈ R 3n A = 1 h Φ i 1 (q (l) ) + D i T 1 n M -1 k(l) , D i T 1 u M -1 k(l) , D i T 1 v M -1 k(l) , 1 h Φ i 2 (q (l) ) + D i T 2 n M -1 k(l) , D i T 2 u M -1 k(l) , D i T 2 v M -1 k(l) , . . . , 1 h Φ i n A (q (l) ) + D i T n A n M -1 k(l) , D i T n A u M -1 k(l) , D i T n A v M -1 k(l) γ ∈ R 3n A = γ i 1 n , γ i 1 u , γ i 1 v , γ i 2 n , γ i 2 u , γ i 2 v , . . . , γ i n A n , γ i n A u , γ i n A v ( 17
)
and the following matrices

D i = D i n , D i u , D i v , i ∈ A(q (l) , ), D = D i 1 , D i 2 , . . . , D i n A , N = D T M -1 D. (18)
Note that the matrix N is positive semidefinite.

In addition, for a vector d ∈ R 3n A , we define by di ∈ R 3 = ( d3 * (i-1)+1 , d3 * (i-1)+2 , d3 * i ). Note that di is a vector, whereas di is a real number component. This convention allows us, after multiplying with M -1 its first equation, to write the problem ( 16) as the conic complementarity problem

(N γ + d) i ∈ -FC i• ⊥ γ i ∈ FC i , i = 1, 2, . . . , n A . ( 19
)
3 Convergence Theory of the Iterative Method

We now describe the structure of projection operators over direct sums of cones. Assume that we have a set of closed convex cones Υ i ⊂ R n i , where the index takes the values i = 1, 2, . . . , n k . We consider the Cartesian product of such cones Υ = n k i=1 Υ i , which we assume is a cone in R n c , that is, that the sum of the dimensions of the element cones satisfies n c = n k i=1 n i . In this section and in the sequel, for a vector x ∈ R n c , we denote by x i , i = 1, 2, . . . , n k its components that satisfy x i ∈ R n i , that is, x = (x 1 , x 2 , . . . , x n k ). Since all the operations we will carry out will be on blocks corresponding to the partition of x into its components x i , there will be no confusion between x i and the components of x. Note that the Cartesian product cone is also a convex cone. Note that we have chosen to use subscripts to denote indices of blocks of the vector x in order to avoid collusion with iteration indices. When particularizing the results to the case of the cone complementarity problem [START_REF] Lotstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF] we will again use superscripts for variables γ pertaining to a contact with index i.

For a convex cone, C ⊂ R m , we denote by Π C (y) the projection of the vector y ∈ R m onto the convex cone C. From the theory of convexity, it follows that the projection has the following properties.

P1 Π C (y 1 ) -Π C (y 2 ) 2 ≤ Π C (y 1 ) -Π C (y 2 ) , y 1 -y 2 , ∀y 1 , y 2 ∈ R m [16][Proposition 3.1.3]. P2 x = Π C (y) ⇔ x ∈ C, y -x ∈ C • , x, y -x = 0 [16][Proposition 3.2.3]. P3 Π Υ (x) = (Π Υ 1 (x 1 ), Π Υ 2 (x 2 ), . . . , Π Υ n k (x n k )) P4 Υ • = n k i=1 Υ i,•
The last two properties are a straightforward application of the properties of convex cones and their projections.

Consider now the symmetric positive semidefinite matrix N . We define the following cone complementarity problem:

(CCP ) s i = (N x + q) i ∈ -Υ i,• , x i ∈ Υ i , x i , s i = 0, i = 1, 2, . . . , n k . ( 20
)
It is immediate that it represents the optimality conditions of the following optimization problem with conic constraints:

(OC) min f (x) = 1 2 x T N x + q T x s.t. x i ∈ Υ i , i= 1, 2, . . . , n k .
The goal of this section is to analyze the following iterative method. We start with an arbitrarily chosen initial point x 0 ∈ Υ . The iterative method is defined by the formula

x r+1 = λΠ Υ x r -ωB r N x r + q + K r x r+1 -x r + (1 -λ) x r , r = 0, 1, 2, . . . , (21)
where λ, ω are parameters that satisfy 0 < λ ≤ 1, ω > 0; for each r, the matrix K r is a strictly block upper triangular or strictly block lower triangular, with blocks corresponding to the partition of the vector x ∈ R n c into the components x i as outlined in the beginning of the section. In addition, B r is a positive diagonal matrix, which is made of identity blocks whose sizes correspond to the same partition of the vector x. We therefore have

B r = ⎛ ⎜ ⎜ ⎜ ⎝ η 1 I n 1 0 • • • 0 0 η 2 I n 2 • • • 0 . . . . . . . . . . . . 0 0 • • • η n k I n n k ⎞ ⎟ ⎟ ⎟ ⎠ , L r = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 0 K 12 K 13 • • • K 1n k 0 0 K 23 • • • K 2n k 0 0 0 • • • K 3n k . . . . . . . . . . . . . . . 0 0 0 • • • 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ , (22)
where

η i > 0, i = 1, 2, . . . , n k , I n i ∈ R n i ×n i , K ij ∈ R n i ×n j , 1 ≤ i < j ≤ n k ,
and we have either that K r = L r , or that K r = L rT . We will use the following assumptions.

A1 The matrix N of the problem (CCP) is symmetric and positive semidefinite. A2 There exists a positive number, α > 0 such that, at any iteration r, r = 0, 1, 2, . . ., we have that B r αI A3 There exists a positive number, β > 0 such that, at any iteration r, r = 0, 1, 2, . . ., we have that (

x r+1 -x r ) T (λωB r ) -1 + K r -N 2 (x r+1 -x r ) ≥ β x r+1 -x r 2 .
To analyze the convergence behavior of the iteration [START_REF] Moreau | Standard inelastic shocks and the dynamics of unilateral constraints[END_REF], we used the same approach as Murty [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF], adapted to the case of general convex cones. We first characterize the solution of the cone complementarity problem in terms of a fixed point of an appropriate mapping.

Theorem 1 Assume that B is a positive definite diagonal matrix with a block structure prescribed in [START_REF] Moreau | Numerical treatment of contact and friction: The contact dynamics method[END_REF]. The vector x ∈ Υ is a solution of the cone complementarity problem (CCP) if and only if it satisfies the following fixed-point relationship:

Π Υ (x -ωB (N x + q)) = x.
Proof By the second property of projections [P2] we have that the vector x satisfies the fixed-point relationship if and only if it satisfies the following relationships:

x ∈ Υ, (x -ωB (N x + q)) -x = s ∈ Υ • ,
x, s = 0.

In turn, these are equivalent to

x i ∈ Υ, (-ωB (N x + q)) i = s i ∈ Υ • , x i , s i = 0, i = 1, 2, . . . , n k .
From the property (P3) of the cones and the fact that the diagonal matrix B has the structure described in [START_REF] Moreau | Numerical treatment of contact and friction: The contact dynamics method[END_REF], such that the blocks corresponding to the components of x are multiples of the identity, it immediately follows that

(-ωB (N x + q)) i = η i ω (-(N x + q) i , ) ,
where we have used the notation from [START_REF] Moreau | Numerical treatment of contact and friction: The contact dynamics method[END_REF]. In turn, this implies that the previously displayed equation is equivalent to

x i ∈ Υ, (-(N x + q)) i = s i ∈ Υ • , x i , s i = 0, i = 1, 2, . . . , n k ,
which is precisely (CCP). The proof is complete.

Theorem 2 Assume that B is a positive definite matrix with the structure described in [START_REF] Moreau | Numerical treatment of contact and friction: The contact dynamics method[END_REF]. Then ∀x ∈ R n c we have that

(Π Υ (x) -x) T B -1 (Π Υ (x) -y) = (Π Υ (x) -x) , B -1 (Π Υ (x) -y) ≤ 0, ∀y ∈ Υ.
Proof From the definition of the total cone Υ we have that Υ = n k i=1 Υ i . Since the matrix B is diagonal with the structure described in [START_REF] Moreau | Numerical treatment of contact and friction: The contact dynamics method[END_REF], we immediately have that

(Π Υ (x) -x) , B -1 (Π Υ (x) -y) = n k i=1 1 η i (Π Υ (x) -x) i , (Π Υ (x) -y) i = n k i=1 1 η i Π Υ i (x i ) -x i , Π Υ i (x i ) -y i .
The last relation follows from the property [P3] of the cones and projections onto them. It is therefore sufficient to show that

x i ∈ R n i ⇒ Π Υ i (x i ) -x i , Π Υ i (x i ) -y i ≤ 0, ∀y i ∈ Υ i . ( 23
)
Using property [P1] of the cones, we have that

Π Υ i (x i ) -y i , Π Υ i (x i ) -y i ≤ Π Υ i (x i ) -y i , x i -y i , ∀y i ∈ Υ i , x i ∈ R n i .
Using the fact that the scalar product is a bilinear form and taking the term from the right to the left with a change sign, we obtain that

Π Υ i (x i ) -y i , Π Υ i (x i ) -x i ≤ 0, ∀y i ∈ Υ i , x i ∈ R n i ,
which proves the equation ( 23) and therefore the theorem. The proof is complete.

Theorem 3 Let {x r : r = 1, 2, . . .} be the sequence of points obtained under the iterative scheme [START_REF] Moreau | Standard inelastic shocks and the dynamics of unilateral constraints[END_REF]. Assume that x 0 ∈ Υ and that the sequences of matrices B r and K r are bounded. Then we have that

f (x r+1 ) -f (x r ) ≤ -β x r+1 -x r 2
for any iteration index r, and any accumulation point of the sequence x r is a solution of (CCP).

Proof The proof is identical to the proof from Murty [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF], where the projection on the positive orthant, +, is replaced with the projection on the cone Π Υ . The only property of the projection that is used is the one from Theorem 2, which holds for the general case as well. We nonetheless include it here for completeness.

Since the initial point satisfies x 0 ∈ Υ and from ( 21), we conclude that x r ∈ Υ , ∀r = 1, 2, . . . From straightforward manipulation it follows that

f x r+1 -f (x r ) = ωB r (N x r + q) T (ωB r ) -1 x r+1 -x r + x r+1 -x r N (x r+1 -x r ) 2 (x r+1 -(1-λ)x r ) λ -x r + ωB r N x r + q + K r x r+1 - (ωB r ) -1 x r+1 -x r + x r+1 -x r N 2 -(λωB r ) -1 -K r x r+1 -x r = λ (x r+1 -(1-λ)x r ) λ -x r -ωB r N x r + q + K r x r+1 -x r T (ωB r ) -1 (x r+1 -(1-λ)x r ) λ -x r + x r+1 -x r N 2 -(λωB r ) -1 -K r x r+1 -
From (21) we know that

x r+1 -(1 -λ) x r λ = Π Υ x r -ωB r N x r + q + K r x r+1 -x r .
We also know that λ > 0. Using these and Theorem 2, we conclude that the first term in the right-hand side of the long equality above is ≤ 0. We therefore have that

f x r+1 -f (x r ) ≤ x r+1 -x r 1 2 N -(λωB r ) -1 -K r x r+1 -x r ≤ -β x r+1 -x r 2 . ( 24
)
The last inequality follows from conditions [A3] and proves the first part of our claim.

Since β > 0, the equation ( 24) implies that f (x r )f x r+1 ≥ 0. Hence {f (x r ) : r = 1, 2, . . .} is a monotone nonincreasing sequence of real numbers. Let x be an accumulation point of the sequence {x r : r = 1, 2, . . .}. Hence, there exists a sequence of positive integers such that the sequence of x r with r belonging to this subsequence of integers converges to x. Since the sequences of B r and K r are bounded sequences of matrices, we can again find a subsequence of the above sequence of positive integers satisfying the property that both the subsequences of B r and K r with r belonging to this subsequence converge to the limits. Let {r t : t = 1, 2, . . .} be this final subsequence of positive integers. Therefore, we can assume that lim t→∞ B r t = B, lim

t→∞ K r t = K, lim t→∞ x r t = x.
In addition, from property (A2) it follows that B is a diagonal matrix with positive diagonal entries. Since f (x) is a continuous function, we have that f (x) = lim t→∞ f (x r t ). Since {f (x r ) : r = 0, 1, . . .} is a nonincreasing sequence of real numbers with a convergent subsequence, f (x r t ), it follows that the entire sequence is itself convergent. This and (24) together imply that 0 = lim t→∞ f (x r t )f x r t +1 ≥ lim t→∞ β x r tx r t +1 2 ≥ 0. From this and the fact that the sequence {x r t } is convergent to x, it follows that

x r t +1 is also convergent to the same limit. These facts imply that

0 = lim t→∞ x 1+r t -x r t = λ Π Υ x r t -ωB r t N x r t + q + K r t ] x r t +1 -x r t -x r t = λ Π Υ x -ω B (N x + q) -x .
So we have that Π Υ xω B (N x + q)x = 0. By Theorem 1, we have that (M x + q, x) is a solution for CCP.

Note that the preceding result does not mean that the sequence will converge, since it is still possible that the sequence will diverge to infinity and have no accumulation point. The proper alternative is related by the following result.

Theorem 4 Under the assumptions of the section, either (a) the sequence x r is bounded, or (b) there exists a 0 = y ∈ Υ, that satisfies N y = 0. In case (a), any two accumulation points z 1 and z 2 satisfy N z 1 = N z 2 .

Proof Assume that case (a) does not hold. Then the sequence x r will have a subsequence x r i that satisfies x r i → ∞, as i → ∞. Consider the sequence

y i = x r i x r i ,
which, being bounded, must have an accumulation point ȳ. We assume, without loss of generality, that the entire sequence y i converges to ȳ ∈ Υ . From the previous theorem, we have that the sequence f (x r i ) is decreasing, and we obtain that

f (x r i ) x r i 2 = 1 2 x r i x r i T N x r i x r i + q T x r i x r i x r i .
Taking the limit as i → ∞, we obtain that ȳT N ȳ ≤ 0. Using assumption [A1], we obtain that ȳT N ȳ = 0. From assumption [A1] it follows that ȳ is a minimum for the function y T N y, and from the optimality conditions it follows that N ȳ = 0. Since ȳ ∈ Υ , case (b) must hold. We have thus proved that the outcome of the iterative method can be only (a) or (b).

Assume now that we are in case (a), and we have two accumulation points z 1 ∈ Υ and z 2 ∈ Υ . From Theorem 3 we have that both z 1 and z 2 are a solution of (CCP), and therefore they satisfy the following relationships.

z 1 ∈ Υ, z 2 ∈ Υ, -(N z 1 + q) ∈ Υ • , -(N z 2 + q) ∈ Υ • .
Since both the cone Υ and the cone Υ • are convex sets, it follows that z 2 + λ z 1z 2 ∈ Υ and that N z 2 +q +λ N z 1 -N z 2 ∈ -Υ • , for any parameter λ ∈ [0, 1]. In turn, from the definition of the polar cone, it follows that

g (λ) = z 2 + λ z 1 -z 2 , Nz 2 + q + λ N z 1 -N z 2 ≥ 0, ∀λ ∈ [0, 1] .
Using the fact that z 2 satisfies z 2 , Nz 2 + q = 0 being a solution of (CCP), we obtain that

λ z 1 -z 2 , Nz 2 + q + λ z 2 , N z 1 -N z 2 + λ 2 z 1 -z 2 , N z 1 -N z 2 ≥ 0, ∀λ ∈ [0, 1]
from which it follows that

z 1 -z 2 , Nz 2 + q + z 2 , N z 1 -N z 2 ≥ 0.
By using the same argument but switching z 1 and z 2 , we obtain that

-z 1 -z 2 , Nz 1 + q -z 1 , N z 1 -N z 2 ≥ 0.
Adding the last two equations, we obtain that

z 1 -z 2 , Nz 2 -N z 1 + z 2 -z 1 , N z 1 -N z 2 ≥ 0,
with which, using the symmetry of the matrix N that follows from assumption A1, we obtain that

2 z 1 -z 2 , N z 1 -N z 2 ≤ 0.
Using again assumption A1, we have that N z 1 = N z 2 , which completes the proof.

Corollary 1 Assume that the friction cone of the configuration is pointed (that is, there does not exist a choice of reaction forces whose net effect is zero). If the relevant parameters satisfy assumptions [A2] and [A3], then the algorithm [START_REF] Moreau | Standard inelastic shocks and the dynamics of unilateral constraints[END_REF] for CCP applied to [START_REF] Lotstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF] produces a bounded sequence, and any accumulation point results in the same velocity solution.

Proof Assume that the sequence x r is produced by the algorithm ( 21) whose parameters satisfy Assumptions [A2] and [A3]. Then, from Theorem 4 there exists 0 = y ∈ Υ such that N y = 0. From [START_REF] Kim | Deep: Dual-space expansion for estimating penetration depth between convex polytopes[END_REF] this implies that Dy = 0.

In turn, from the definition of D in ( 18) and Subsection 2.3, this implies that there exist nonzero constraint feasible impulses that produce a zero net effect on the system. This contradicts the assumption that the friction cone is pointed [START_REF] Stewart | An implicit time-stepping scheme for rigidbody dynamics with inelastic collisions and Coulomb friction[END_REF]. Therefore boundedness of the iteration sequence and the existence of an accumulation point are assured. Uniqueness of the velocity follows from the second part of Theorem 4, since N z 1 = N z 2 and the definition of D in [START_REF] Kim | Deep: Dual-space expansion for estimating penetration depth between convex polytopes[END_REF] implies that Dz 1 = Dz 2 , which, in turn, from [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms[END_REF] implies that the velocity solution is unique. The proof is complete.

Note that N is symmetric positive semi-definite and therefore assumption [A1] is satisfied. Assumption [A2] is easily satisfied, whereas assumption [A3] can be satisfied by a trial-and-error approach, whereas if the iterates x r+1 and x r do not satisfy [A3], than the parameter ω is decreased by factor of two, and the iterate x (r+1) is recomputed. It is immediate that such a strategy can decrease the parameter ω only a finite number of times.

Implementation

The CCP method proposed here can be applied to the simulation of multibody systems with a large number of parts and contacts because, where an upper limit on the number of iteration is enforced, the iteration ( 21) can run in O(n) space and O(n) time.

Previous sections showed that generic multibody problems with frictional contacts, expressed with the system ( 5)-( 8), embed the cone complementarity problem [START_REF] Lotstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF]. Hence, the iterative method ( 21) can be used to solve such convex CCP, because ( 19) is equivalent to the more general problem [START_REF] Marques | Differential Inclusions in Nonsmooth Mechanical Problems: Shocks and Dry Friction, volume 9 of Progress in Nonlinear Differential Equations and Their Applications[END_REF] where one considers the specific case of three-dimensional cones Υ i . That is, for the ith friction cone Υ i we have that n i = 3 and that there is an associated vector with a normal and two tangential reactions:

γ i = {γ i n , γ i u , γ i v }.
The complete vector of unknown scalar reactions is γ ∈ R 3n A . From Section 3, we have that n k = n A and n c = 3n A .

Given [START_REF] Jourdan | A Gauss Seidel like algorithm to solve frictional contract problems[END_REF], the final time-stepping scheme can be seen as a sequence of three main operations: a CCP problem that finds unknown reactions γ (25a), a linear application (25b) that gives the new speeds v (l+1) , and a position update (25c):

(N γ + d) ∈ -Υ o ⊥ γ ∈ Υ (25a) v (l+1) = M -1 k + Dγ (25b) q (l+1) = q (l) + hΓ (q (l) )v (l+1) . (25c)
The biggest computational overhead is caused by the first problem, that is, the CCP (25a). In fact, (25c) is immediate, and (25b) can be computed quickly because in most cases the matrix M is diagonal and its inverse M -1 can be precomputed easily.

We recall that N = D T M -1 D. The full D matrix can be partitioned in n A vertical blocks D i ∈ R m×3 , each pertaining to the corresponding ith cone. We also recall that, from [START_REF] Kim | Deep: Dual-space expansion for estimating penetration depth between convex polytopes[END_REF], we have that

D = D 1 |D 2 |...|D n A .
Using [START_REF] Jourdan | A Gauss Seidel like algorithm to solve frictional contract problems[END_REF], we can rewrite the term r from [START_REF] Jourdan | A Gauss Seidel like algorithm to solve frictional contract problems[END_REF] in a more compact form:

d = D T M -1 k + b. (26)
The convergence theory about the iterative scheme ( 21) leaves some degrees of freedom in choosing η i values that build the diagonal blocks of the iteration matrix B. A trivial choice could be to use the same η i = ξ value for all diagonal blocks, that is, B = ξI, and then use the overrelaxation parameter ω to control the convergence. However, setting the same value for all η i may slow convergence in systems with large mass ratios, even with an optimal ω. A more practical approach, which copes better with systems affected by uneven masses, is to use η i = 1/ḡ i , where ḡi is the average of the diagonal values of the ith block of the N matrix. We note that ḡi can be computed easily from the trace of the 3

× 3 matrix D i,T M -1 D i , as ḡi = Trace(D i,T M -1 D i ) 3 . (27)
Also the K matrix in ( 21) can be chosen freely, within the convergence limits posed by assumptions [A1]-[A3]. Among the most noticeable options, we note the case where K = 0, which results in a scheme like a projected Jacobi, or the case where K is built by using the lower blocks of N , so that K i,j = D i,T M -1 D j , where 1 ≤ j < i ≤ n A . In practical terms this means that, as soon as computed, a triplet γ i with three reaction values will be used also for computing the following γ i+1 triplet, and so on for all i, without needing to finish a single iteration, which results in a Gauss-Seidel-type iteration. Numerical tests show that this last option, similar to a projected SOR scheme with immediate update of unknown vector, converges faster than the case of K = 0. Hereafter, we will assume that such a kind of K matrix is used. Another choice, that we do not explore here, is the one of having the matrix K be block-diagonal with block lower triangular blocks, which would be equivalent with block Jacobi, where for each block we do a Gauss-Seidel-type iteration. The latter is suitable for a parallel iteration with low communication overhead.

We recall that the matrix N is a product of large matrices; N = D T M -1 D, and it is full even if D and M are sparse. For systems with a large number of contacts, the size of N would be prohibitive and clearly would not satisfy the goal of O(n) space complexity. To this end, direct multiplication of vectors and matrices in (21) must be avoided; otherwise the effort and the space requirement would be superlinear in the number of constraint.

For the reasons above, a scheme that does not need N computed explicitly has been developed, exploiting the sparsity of M and D. Also the K matrix does not need to be explicitly built, if we adopt the above mentioned choice of K as the upper block-structure of N . These considerations lead to the following implementation of the rth step of the iteration [START_REF] Moreau | Standard inelastic shocks and the dynamics of unilateral constraints[END_REF], expressed as an inner loop with index i = 1 . . . n A on all n A friction cones Υ i :

δ ir = γ ir -ωη i D i,T M -1 i-1 z=1 D z γ z,r+1 + + n A z=i D z γ zr + ki + b i (28)
γ i,r+1 = λΠ Υ i δ ir + (1 -λ)γ ir . (29)
In the case of friction in three-dimensional space, the implementation of the projection operator

Π Υ i (δ i ) : R 3 → R 3 is straightforward.
For improved performance, some operations can be computed at the beginning of the iteration because their outcome would remain unchanged. In detail, we introduce the m × 3 matrix s i = M -1 D i and the 3 × 3 matrix

g i = D i,T M -1 D i .
Considering the optimizations above, we can express the final CCP algorithm with the following pseudocode:

Algorithm 1 1. For i = 1, 2, . . . , n A compute the m × 3 matrices s i = M -1 D i and 3 × 3 matrices g i = D i,T s i . 2. For i = 1, 2, . . . , n A , compute η i = 3/Trace(g i ).
3. If warm starting with some initial guess γ * , initialize reactions as γ 0 = γ * , otherwise γ 0 = 0.

Initialize speeds

: v = n A i=1 s i γ 0 + M -1 k. 5. For i = 1, 2, . . . n A , perform the updates δ i,r = γ i,r -ωη i D i,T v r + b i ; γ i,r+1 = λΠ Υ δ i,r + (1 -λ)γ i,r ; Δγ i,r+1 = γ i,r+1 -γ i,r ; v := v + s i T Δγ i,r+1
. 6. Repeat the step 5 by looping on the list of contacts in backward direction, if symmetric updates are desired. 7. r := r + 1. Repeat from 5 until convergence, or until r > r max .

The iterations, usually stopped when an approximation threshold has been reached, can be also prematurely aborted when r exceeds a limit r max on the maximum number of iterations if the simulation must meet hard-realtime requirements.

With minimal changes to the Π Υ (•) operator, the proposed method can easily adapted to the case of friction in 2D or the case of generic unilateral constraints. Also, without major modifications to the main scheme, classical bilateral constraints can be added. In that case, the cones Υ i can be taken to be R l i , where l i is the dimension of the force vector associated to a joint and the projection on such a cone is the identity. Therefore the scheme is modified by simply replacing [START_REF] Pang | Differential variational inequalities[END_REF] with

γ i,r+1 = λδ ir + (1 -λ)γ ir . ( 30
)
In our simulations, we chose ω = 1 and λ = 1, except for the test of convergence of the residual with respect to ω. We cannot guarantee a priori that this will satisfy condition [A3], but it did for all our simulations. In addition, the matrix sequences K r and B r were constant. We can therefore claim that Theorem 3 does apply and, since the sequence did not diverge (and was in fact convergent), any accumulation point is a solution of the cone complementarity problem (25a). In addition, Theorem 4 is applicable to show that any accumulation point has the same velocity solution. It is difficult to verify numerically the condition of Corollary 1 (the pointed friction cone assumption). Nonetheless, boundedness of the iterates was observed in all cases. In addition, our proofs of the theoretical results allow for similar conclusions if ω varies from iteration to iteration. Therefore, we could ensure that at some iteration the appropriate ω is chosen after decreasing its value a few times until assumption [A3] holds. It can be shown that if the value of ω is halved each time [A3] does not hold and the respective iteration is rejected, then [A3] will eventually be satisfied after a finite number of steps. In our experiments, however, the values we have chosen for ω and λ have worked for all iterations without need of further adjustment.

Finally, the overall scheme [START_REF] Pang | Convergence of time-stepping method for initial and boundary-value frictional compliant contact problems[END_REF] for advancing a single time step can be expressed with the following pseudocode: Algorithm 2 1. Set t = 0, step counter l = 0, provide initial values for v (l) and q (l) . 2. Perform collision detection between shapes of bodies, obtaining n A possible contact points within an distance. For each contact point, compute D n , D u , D v and residual Φ(q). 3. For each body, compute forces f t (q, v, t), then compute k with [START_REF] Haug | Dynamic mechanical systems with coulomb friction, stiction, impact and constraint addition-deletion[END_REF] and b with (17). 4. Use Algorithm 1 to obtain unknown impulses γ and speeds v (l+1) for the CCP problem. 5. Update positions using q (l+1) = q (l) + hΓ (q (l) )v (l+1) . 6. Increment t := t + h, l := l + 1, and repeat from step 2 until t > t end

We remark that choosing a proper value for the collision envelope is not trivial. If a very small or zero value is used, contacts will enter the CCP solver only when it is too late and some amount of interpenetration will be unavoidable: this will affect negatively the stability of the method. On the other hand, if too large values are used, the collision detection algorithm will return too many potential contacts which will waste computational resources for the CCP solution and which could occasionally create troubles with convex shapes: this will decrease the efficiency and the robustness of the method. We experienced that a simple yet efficient heuristic is to choose as the maximum distance which can be spanned by whatever point of the moving body, given its speed.

Examples

We present the results of our algorithm on two granular materials applications. For the larger simulation, the number of impulse variables exceeded 400,000.

Size segregation in a shaker

The first example is meant as a benchmark to evaluate the performance of the solver when dealing with many contacts with friction. A rectangular box is filled with spheres; then the box is shaken by means of an articulated suspension and a crank mechanism (Fig. 2). When large objects are mixed with the spheres, a phenomenon called vibration-induced size segregation moves larger objects on top: this effect can be observed also in our simulations.

Different parameters have been tested, for example repeating the simulation with a varying number of spheres up to 1500. In all cases the mass of the spheres is m = 0.01kg, their diameter is d = 26mm, the friction coefficient is μ = 0.3, and the time step is h = 2π/50Ωs, with Ω rad/s being the frequency of the crank. The amplitude of the vibration has been tested up to A = 10mm.

Plotting of ||Δγ|| (Fig. 3) during the iteration of the algorithm shows the convergence of the method for varying values of the overrelaxation factor ω.

Figure 4 shows how the CPU spends time in various parts of the simulation algorithm. For this benchmark, a shaker with 1,000 rigid bodies was simulated with an upper limit of 40 iterations for the CCP solver. One can easily see that the solution of the CCP is the bottleneck in the entire simulation process, while collision detection and other tasks (time integration, Jacobian update, etc.) are less CPU-intensive. In this example, the 40-iteration limit was enough to keep the feasibility errors at negligible levels (max. interpenetration || P n || < 0.002 d). However, if lower precision is acceptable as in case of virtual reality or real-time applications, fewer iterations can be used, thus reducing the CCP timings to levels which are comparable to the collision detection timings.

To show how the number of iterations affect the precision of the solution to the complementarity problem, in Fig. 5 we report the maximum error in terms of speed violation || V n || in contact constraints, during 300 time steps of simulation. Speed is measured in d/s, where d is the diameter of the spheres. Similarly, we report in Fig. 6 the maximum position error || P n || in contact constraints, that is, the maximum interpenetration. The error is measured in d units. One can see that, despite the large number of objects in contact, acceptable precision can be obtained also with a moderate number of iterations.

By performing a set of six shaker simulations with an increasing number of objects, hence for increasing numbers of contacts, one can obtain a graph as in Fig. 8, which shows how the CPU effort grows linearly with the number of frictional contacts. Here, for each simulation, CCP timings have been recorded after ten seconds of transient, when spheres are at steady state and form a dense packing, because this is a nontrivial configuration that requires significant CPU efforts. The linear-time complexity is a consequence of the loop in the fifth step of the algorithm, which is O(p) with p reaction forces γ i if a maximum number of iterations is enforced (40 iterations in this example). Note that the fourth step of the algorithm, performing a computation that is linear in terms of number of rigid bodies, has a moderate or negligible impact on overall performance. Despite the fact that the theoretical complexity of the algorithm for fixed number of iterations is linear, some deviation from linearity can be experienced in complex applications when large amounts of contacts are simulated, because CPU cache misses can become more frequent as the memory access starts to become more and more intense.

Figure 7 show that only a portion of the potential contacts will be active (i.e., with nonzero reaction force) after the CCP solution. Since the computational effort is proportional to the number of potential contacts entering the CCP solver, regardless of the active/inactive outcome, a proper collision detection algorithm should take care to report the smallest number of potential contacts, that is, only the surface pairs that may give interpenetration in a single time step given the actual state of bodies. This precaution would keep the active contacts as a relatively fixed percentage of the number of potential constraints.

Granular flow from a silo

The numerical method proposed in this article can be used to simulate dense granular flows in silos. This problem arises in many engineering applications, most noticeably in the development of the promising fourth-generation uranium-based, graphite-moderated, helium-cooled very high temperature Fig. 9 Frames from the simulation of 36,000 rigid bodies with frictional contacts flowing from a three-dimensional funnel. On a T2600 2 GHz processor, each solution of the CCP problem (nearly half a million of variables, including constraint multipliers and speeds) with 140 iterations, took 19 s of CPU time on average. nuclear reactor, where thousands of graphite fuel pebbles drain very slowly in a continuous refueling process [START_REF] Gougar | Advanced core design and fuel management for pebble-bed reactors[END_REF].

Pebble flow in such pebble-bed reactors (PBRs) is not easily accessible to experiments, and no reliable continuum model is yet available for analytical approaches. These facts motivate the development of fast numerical methods. Simulations of PBR reactors have been recently performed with DEM discrete-element methods [START_REF] Rycroft | Analysis of granular flow in a pebble-bed nuclear reactor[END_REF]; however, the DEM approach is based on a stiff spring-dashpot contact model which requires a very small time step in order to guarantee the stability of the integration. Conversely, the method proposed here can enforce rigid contacts without the need of artificial stiffness; hence larger timesteps can be used. For example, the flow simulation of Fig. 9, representing 11 seconds of drainage from a silo 3.5m wide with 36,000 uranium-graphite spheres with d = 0.06m and friction coefficient μ = 0.6, exploited a timestep h = 0.01s that is three orders of magnitude larger than the timestep required by the DEM method in [START_REF] Rycroft | Analysis of granular flow in a pebble-bed nuclear reactor[END_REF]. The simulation took about four hours to complete for 5 seconds of simulated time, with a penetration error comparable to the one in the size segregation case. But timing is perhaps less relevant since it depends on items such as cache management that can vastly change with different optimization than the fact that the simulation completed with low penetration error for a fixed (and relatively small) number of iterations, 140, for a very high density configuration. The maximum number of contacts for which the problem was benchmarked was almost half a million, which in turn resulted in more than two millions of variables for the CCP. This is a promising approach to the simulation of full-scale reactors and other large granular flow problems, though further tests are needed to determine whether the maximum penetration error does not increase with an increasing number of uranium-graphite spheres.

Conclusions

Aiming at a linear-time solution of dynamical systems with thousands of constraints and contacts, we have presented a novel method for solving the cone constrained subproblems that appear in a time-stepping approach recently proposed in [START_REF] Anitescu | Optimization-based simulation of nonsmooth rigid multibody dynamics[END_REF]. The method has the flavor of a Gauss-Seidel with overrelaxation and is proven to converge under fairly standard assumptions about the configuration of the system.

We implemented this method into the HyperOctant library of our multibody project, Chrono::Engine [START_REF] Tasora | Chrono::engine project[END_REF]. Our method is able to handle large simulations with tens of thousands of colliding rigid bodies and hundreds of thousands of constraint impulse variables and scales well in this range. In previous work [START_REF] Anitescu | A fixed-point iteration approach for multibody dynamics with contact and friction[END_REF][START_REF] Tasora | Un nuovo metodo del simplesso per il problema di complementarit lineare mista in sistemi multibody con vincoli unilateri[END_REF] we have shown that simplex-like methods do not scale well for systems in configurations of the type solved here. In future work, and as appropriate software packages become available, we will carry out comparisons with interior-point methods for optimization problems with conic constraints.

Because of the low computational overhead of our method, we foresee that it could be endorsed even in the emerging application fields of physical engines for videogames and virtual interactive environments, which can exploit the benefits of the method for real-time performance.

Fig. 1

 1 Fig. 1 Examples of multibody systems with many unilateral contacts.

Fig. 2

 2 Fig.2Shaker benchmark, for a 1500-sphere case, showing vibration-induced granular segregation of large objects. After ten seconds of simulation three black intruder particles, originally placed at the bottom of the shaker, will rise at the top surface.

2 Fig. 3 Fig. 4

 234 Fig.3Convergence of Δγ r for varying ω, for a sample time step in the 300sphere benchmark.

10 Fig. 5 10 Fig. 6 Fig. 7

 1051067 Fig.5Maximum speed violation in constraints, for the 300-sphere benchmark.

Fig. 8

 8 Fig. 8 Average CPU time used to compute a step of simulation, as a function of the number of contacts.

Though the original Coulomb model distinguishes between static μ s and kinetic μ k friction coefficients, where usually the kinetic coefficient is slightly lower than its static counterpart, in this work we consider both to have the same value μ because the difference is not relevant for the discussion and suffices to say that a proper algorithm might adjust the friction coefficient adaptively during the simulation, depending on the slipping speed, to match complex nonlinearities in μ as a function of speed.

Acknowledgments

We thank Paul Tseng for technical discussions concerning block coordinate descent methods. We are profoundly grateful to the two referees and to Dan Negrut for their constructive comments. Mihai Anitescu was supported by Contract No. DE-AC02-06CH11357 of the U.S. Department of Energy.

A Notations

A.1 Multibody system n: number of bodies.

m = 6n: dimension of the position state vector.

-M : Mass matrix, positive definite, of size m × m.

q: vector of generalized positions of dimension m.

v: vector of generalized velocites of dimension m.

t: time of the system.

h: time step used by the time-stepping scheme.

f t (t, q, v), f e (t, q, v), f c (q, v): the total, external, and, respectively, Coriolis forces acting on the system, vectors of dimension m.

p: number of contact constraints that can become active (no more than

).

-A(q l , ) set of -active contact constraints, a set with no more than p elements.

n A : dimension of set A(q l , ).

-D: aggregate matrix of normal and tangential directions at the contact in generalized coordinates, a matrix of dimension m × 3n A .

A.2 Contact and friction model

-Φ(q): the gap function, which indicates whether a contact constraint is active.

-n, t 1 , t 2 : the normal and tangential vectors at a contact, three-dimensional vectors. -F N , F T : the normal and tangential force at a contact, three-dimensional vectors.

-v T : the tangential velocity, a three-dimensional vector. -FC: the Coulomb friction cone, a subset of a three-dimensional space.

A.3 Cone Complementarity Problems

n k : number of cones Υ i whose direct sum give the total constraint cone of the cone complementarity problem.

n i : dimension of the vector space R n i in which the component Υ i is embedded.

n c = P n k i=1 n i : dimension of the unknown vector of the cone complementarity problem.

-Υ : the total constraint cone of the cone complementarity problem, a subset of -B: a block diagonal matrix, with multiple of identity blocks.

ω, λ: scalar parameters of the iterative scheme.

The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory ("Argonne"). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government.


John Rieffel
email: rieffelj@union.edu

Jean-Baptiste Mouret
email: jean-baptiste.mouret@inria.frnicolas.bredeche@isir.upmc.fre.haasdijk@vu.nl

Nicolas Bredeche

Evert Haasdijk

Introduction to the Evolution of Physical Systems

culmination of a series of workshops organized around the topic

A Brief Survey of the Evolution of Physical Systems

The Evolution of Physical Systems has its roots in the embodiment philosophy of Rodney Brooks, who famously said "the world is its own best model" [START_REF] Rodney | Elephants don't play chess[END_REF]. Brooks's emphatic critique of the Symbol System Hypothesis argued that relatively simple systems, developed and grounded in suitably complex real-world environments, can lead to the emergence of complex behaviors. Contemporary approaches to the Evolution of Physical Systems aim at automating Brooks' vision by exploiting evolutionary algorithms on systems that "live" and "evolve" in the real world, that is, whose behavior (and possibly their form) are as grounded as possible in the environment.

Embodied evolutionary algorithms first gained prominence with the "Sussex Approach" of Harvey et al. [START_REF] Harvey | Evolutionary robotics: the sussex approach[END_REF] in the 1990s. The Sussex gantry-based robot [START_REF] Cliff | Evolving visually guided robots[END_REF][START_REF] Harvey | Seeing the light: artificial evolution, real vision[END_REF] driven by a neural network was capable of robust behaviors in a noisy real-world environment. Their motivation for using the real-world rather than a simulated environment stemmed largely from the fact that in that era of limited computational power, it was faster and more efficient to test networks in situ than it was to develop and run a simulator that realistically modeled sensor noise. This notion that it is sometimes more practical (and computationally efficient) to dispense with simulation entirely, and evolve behaviors directly in the real world, persists at the heart of the Evolution of Physical Systems.

Two other results of the Sussex group of particularly noteworthy. First is Adrian Thompson's "Silicon Evolution" of an FPGA [START_REF] Thompson | An evolved circuit, intrinsic in silicon, entwined with physics[END_REF][START_REF] Thompson | Silicon evolution[END_REF]. The chosen task was to evolve a circuit capable of discriminating between 1kHz and 10kHz signals without the use of a clock. By evolving in silicon, rather than in simulation, Thompson was able to generate nearly inscrutable solutions that exploited the analog (and non-simulatable) nature of FPGAs. This idea of using embodiment to arrive at "novel surprise" solutions is a further guiding principle of the Evolution of Physical Systems. The second is Jakobi's essential work on the "Reality Gap" [START_REF] Jakobi | Noise and the reality gap: The use of simulation in evolutionary robotics[END_REF], which highlighted the difficulty in transferring results evolved in simulation into the real world. Jakobi notes that the best results emerge when "the noise levels of the simulation have similar amplitudes to those observed in reality", and points out that these approaches become less feasible as environments and sensors become more complex.

The practice of evolving outside of simulators grew from there. Floreano and Mondada [START_REF] Floreano | Automatic creation of an autonomous agent: Genetic evolution of a neural network driven robot[END_REF] were able to evolve robust wall-following behavior in a tethered Khepera robot through embodiment.

Watson and Ficici [START_REF] Richard A Watson | Embodied evolution: Distributing an evolutionary algorithm in a population of robots[END_REF] took embodiment a step further by simultaneously embedding an entire population of evolving "Tupperbots" into a shared arena -and demonstrated how the behaviors that emerged through this process were qualitatively different, and more effective, than those evolved in a simulated environment (their use of an electrified floor as an alternative to tethers was particularly innovative). Zykov et al. evolved dynamic open-loop gaits on a large hexapod robot, remarking upon the amount of labor required to reset the large robot between physical trials. [START_REF] Zykov | Evolving dynamic gaits on a physical robot[END_REF] Several of our own research areas have been in this realm of problems which are easier to test physically than they are to realistically simulate. In their work on "Evolutionary Fabrication" [START_REF] Rieffel | Evofab: a fully embodied evolutionary fabricator[END_REF][START_REF] Kuehn | Automatically designing and printing 3-d objects with evofab 0[END_REF] Rieffel and colleagues ran an embodied evolutionary algorithm on a 3D printer, learning to produce a variety of shapes out of viscous extruded materials such as silicone. These efforts led to novel shapes (and novel means of producing the shapes) that no human would have designed.

Similarly, and more recently, Rieffel's work on the evolution of vibration-based gaits for tensegrity robots [START_REF] Rieffel | Soft tensegrity robots[END_REF] relied purely on physical trials, and led not only to fast and dynamically complex gaits, but also a state-machine based controller for a tensegrity robot able to chase targets [START_REF] Khazanov | Exploiting dynamical complexity in a physical tensegrity robot to achieve locomotion[END_REF].

While not using a simulator at all has many benefits, evolutionary robotics can take advantage of both simulators and real experiments to accelerate the evolutionary process while keeping a direct link with the real world. A first idea is to start evolution in simulation then finish it on the real robot(s); this is, for instance, what Lipson and Pollack did to evolve the morphology and controllers of 3D-printed robots [START_REF] Lipson | Automatic design and manufacture of robotic lifeforms[END_REF]. However, this approach means that designs for the real world are variants of those evolved in simulation, that is, robots are more grounded in the simulation than in the real world. A second, classic idea is to use data collected from tests on the real robot to improve the simulator, ideally to make the reality gap so small that it becomes unnoticeable [2,[START_REF] Cristóbal | Combining simulation and reality in evolutionary robotics[END_REF]. Nevertheless, to do so, simulators needs to generic enough to model almost anything that could happen, which, in turns, requires acquiring a considerable amount of data in the real world.

More recently, a weaker version of this idea has emerged: instead of trying to correct the simulator, one could use a few tests in the real world to learn to avoid the behaviors that exploit phenomena that are inaccurately simulated. Following this concept, Koos et al. managed to evolve controllers for 6-legged robots and Kheopera like-robot with only a dozen of tests in the real world [START_REF] Sylvain Koos | Fast Damage Recovery in Robotics with the T-Resilience Algorithm[END_REF][START_REF] Sylvain Koos | The Transferability Approach: Crossing the Reality Gap in Evolutionary Robotics[END_REF][START_REF] Mouret | Crossing the reality gap: a short introduction to the transferability approach[END_REF].

Despite numerous advances in both robotics and technology, Rodney Brooks' prescient words about the crucial role of embodiment in robotics remain as true as ever, and the Evolution of Physical Systems is as exciting as ever.

The Papers

We the editors are delighted to present the papers that were selected to be published in this issue.

Collectively they provide an inspiring snapshot of the state of the Evolution of Physical Systems.

Fully Autonomous Evolution of Physical Machines

Jacobsen and colleagues describe an experimental setup that removes the burden of manually setting the environment (e.g.: moving the robot(s) and object(s) back to their initial positions and orientations) from the human supervisor. This is particularly relevant in evolutionary robotics, where a large number of trials are generally required, all starting from similar initial conditions.

The authors present a proof-of-concept, tested with a real robotics setup in an arena of 1m per 0.5m, which makes use of an industrial robot arm, a camera and fiducial markers that can be sticked on top of robots and objects alike. These fiducial markers are used both for monitoring the experiments, by tracking positions and orientations, and as visual anchors for the robot arm to grab and move back robots and objects to their initial positions. Beyond evolutionary robotics, the proposed setup may also be used for other robotics setups involving multiple trials (e.g.: reinforcement learning, embodied evolution, ...), would it be for resetting the experimental setup, as demonstrated in this paper, or for other more dynamical interventions (e.g.: by conditionally moving objects around or transporting robots during the course of an experiment). To ensure their system can be adopted and easily deployed by others, the authors have included the source code and a lengthy description of the technical aspects of their setup.

Abstraction, Sensory-Motor Coordination and the Reality Gap in Evolutionary Robotics

In their contribution, Scheper and de Croon describe how choosing an appropriate level of abstraction can improve the transfer of evolved solutions across the reality gap. The task domain they have have chosen is a decentralized formation flying task by a group of three small airborne vehicles controlled by neural networks. They compare the performance of a low-level controller that has precise control over vehicle dynamics (pitch rate, pitch angle, etc) against a more high level controller that uses a more abstract velocity set-point. After evolution (in simulation) the more abstract-level controller was better able to survive transfer into the real world. Scheper and de Croon suggest that the higher-level controller is better able to perform a form of environmental exploration by inherently modeling the behavior of the other vehicles in the swarm. As they point out, an advantage of more abstracted controllers is that it reduces the necessity of high-fidelity simulation -an important advantage for evolutionary robotics tasks. We were impressed by the insights of this paper in regards to the trade-offs between simulator fidelity and real-world behavior.

The flight videos on YouTube are also quite impressive.

Evolutionary Developmental Robotics: Improving Morphology and Control of Physical Robots

The work of Vujovic and colleagues is very much in the spirit of this special issue. They introduce a novel robot-building robot composed of a hot-melt extruder and an articulated arm, that acts as a sort of 3D printer, capable of fabricating the legs of simple modular robots. An evolutionarydevelopmental approach is used to "grow" the designs of the legs, and once printed, the locomotive abilities of the printed robot are automatically evaluated in situ. This type of physically embodied integrated design-fabricate-test loop, in which one robot evolves and assembles other robots, is very much in the vein of Lipson and Pollack's GOLEM robots [START_REF] Lipson | Automatic design and manufacture of robotic lifeforms[END_REF], moving the field one step closer to 3D printed robots that can walk themselves off the printer.

On Design Mining: Coevolution and Surrogate Models

Preen and Bull present an interesting theoretical analysis is followed up with an actual physical case study. Evaluation of physical artifacts is commonly a time-consuming aspect of the Evolution of Physical Systems, so a thorough analysis of methods to minimize the number of required evaluations can greatly boost applicability. Preen and Bull employ the abstract NKCS model of coevolution to achieve just this: it enables informed sampling of candidate solutions in a co-evolutionary scenario where various interacting components of a physical system in combination with a neural net-based surrogate model. A successful case study to develop a heterogeneous set of vertical wind turbines for wind farms shows the applicability and relevance of their method to evolve designs for physical artifacts without relying on simulation.

Real-World Evolution of Robot Morphologies: A Proof of Concept

For obvious practical reasons, most of the work in the Evolution of Physical Systems has been focused on the evolution of controllers for existing robots, that is, for robots with a fixed morphology. Nevertheless, most of those who are interested in the field dream of a much more ambitious vision: robots whose morphology and brains can evolve on their own, much in the same ways as animals evolve.

In their contribution, to appear in a forthcoming issue, Jelisavcic et al. investigate how such a vision could be achieved with a population of robots and 3D printing technology. Using the RoboGen system of 3D-printed modules [START_REF] Auerbach | Robogen: Robot generation through artificial evolution[END_REF], they first designed and 3D-printed 2 robots. In a second step, they put these two robots in a arena so that they can exchange their genome to create the genotype of a third robot, which they also 3D-printed. This first step of an evolutionary experiment is an interesting stepping stone to identify the challenges of setting up experiments in which a population of robots would evolve autonomously.

Looking Forward

Overall, these papers provide valuable contributions to the field of the Evolution of Physical Systems, demonstrating that many of the current challenges of the field are more technological than scientific: to design machines that evolve in the physical world, we first need machines that can self-reproduce or, at least, autonomous machines that can in turn help other machines to evolve.

And until high fidelity computationally efficient simulators are developed, the role of physical embodiment in Evolutionary Robotics remains essential.


Guest Editorial of the Special Issue: Manufacturing Plant Control Challenges and Issues

We are pleased to introduce this special issue of the IFAC journal on Control Engineering Practice which contains improved selected papers of the 11 th IFAC INCOM'04 Symposium on Information Control Problems in Manufacturing held in Brazil.

The main purpose of INCOM'04 has been to point up international research and developments dealing with all the applications of automation, information and communication technologies in order to control and to manage the manufacturing plant within the e-enterprise. This general scope involves all methodological and technological aspects to digitally control with more agility the entire manufacturing chain, from supply and design through manufacturing, to maintenance and service, over the whole product and processes life cycle.

The introductory paper by (Morel, Valckenaers, Faure, Pereira, Diedrich) presents current key problems/applications and the recent major accomplishments/trends addressed by the three working groups of the IFAC TC for Manufacturing Plant Control issues as well as by the main sponsor of two workshops on Intelligent Manufacturing Systems and Intelligent Assembly & Disassembly and a Symposium on Information Control Problems in Manufacturing.

Then, (Neumann) presents "what is going on communication in industrial automation" discussing how communication networks and the Internet are decisively impacting the manufacturing automation domain. This general problem of automation over and with networks leading to a distributed automation system is handled by the working group on Networked Controlled Manufacturing Systems which has contributed to making TC 5.1 proactive since 2001 in this important application field by investigating the following aspects:

. Design, programming, operation and diagnosis of automation behaviour in a distributed environment, . Industrial informatics methods for automation device and system development, . System integration models and technologies for communication connected devices (e.g. for configuration and parameterisation), . Use of heterogeneous (industrial and private) networks, so called Virtual Automation Networks (VAN), to provide automation relevant Quality of Services, . Life cycle aspect for distributed automation systems.

A first rational issue of this networked automation is to put into question the Hierarchical/Integrated vision of Enterprise-wide control for a more Interoperable/Intelligent one by postulating that the customized product should be the 'controller' of the manufacturing enterprise resources. This is made possible by the related explosion of enabling information technologies among which we can cite infotronics technology such as RFID. (Parkilad and McFarlane) investigate the role of this product information in end-of-life decision making in order to ensure the coherence between the physical and information flows all along the product life cycle. These novel solutions for manufacturing control systems that merge know-how from control engineering, software engineering and complex systems/artificial life research are investigated by the working group on Intelligent Manufacturing Systems. This results in novel designs promising scale-ability, reusability, integrate-ability and robustness. Recent efforts focus on the validation of these concepts through realistic testing campaigns on industrial test cases as stated by (Marik and Lazansky) for industrial applications of agent technology. To this end and in anticipation of the availability of a generic service as investigated by the European IMS Network of Excellence, simulation-based benchmarking of production control schemes for complex manufacturing systems addressed by (Mönch) remains the mean to prove the efficiency of solutions before their deployment for practical purposes. On a more fundamental level, the concepts of emergent and self-organizing systems, often in designs inspired by biological systems, are posing new engineering perspectives, such as the engineering perspective on the supply network control problem addressed by (Saint Germain, Valckenaers, Verstraete, Hadeli, Van Brussel) by extending the HMS paradigm for inter and intraenterprise logistics issues. On an organizational level, a continuing collaboration of this WG with the IFIP WG 5.7, the European IMS Network of Excellence and the Brazilian MANET Network of Research contributes to making TC 5.1 proactive to unite the scientific IMS community.

Another rational issue addressed by (Johnson) is to test the role of formal methods for improving automation software dependability in order to face the tolerance to faults of these networked and software-based control systems. As examples of formally checking their behaviour, (Flordal, Fabian, Akesson, Spensieri) deal with the automatic model generation and PLC-code implementation for coordination of industrial robot cells and (Kwiatkowska,Norman, Parker) with their dependability analysis by probabilistic models. This general problem of the dependability of manufacturing systems driven and supervised by hardware and software controllers that most of the time are communicating through networks is addressed by the working group on Dependability of Manufacturing Systems which has contributed to making TC 5.1 proactive since 2002 by investigating the following aspects:

. Formal or formalized description techniques for design, implementation and validation of system components and communication systems, . Manufacturing systems formal design techniques, . Safe and secure middle ware based on COTS components, . Dependability in networked automation . Safety and security algorithms, . Modes management.

Finally, this advanced manufacturing plant-wide automation is posing new education and training challenges to develop system thinking in order to lead the stepwise abstraction and concretization of technical system complexity towards the emerging ambient intelligence world. (Bruns and Erbe) propose a learning environment merging on-site and remote components into a cooperative systems engineering process in order to bridge reality and virtuality over only idealized computer-simulation learning situations.

We would like to thank the authors for their contributions and the referees for the time they put in reviewing all the papers. Thanks to Professor Ian Craig, Editor-in-chef of the IFAC journal on Control Engineering practice for his strong support.

Pr. Gérard Morel Chair of International Program Committee

Pr. Carlos Eduardo Pereira General Scientific Chair