Olivier Goury

Christian Duriez

Fast, generic and reliable control and simulation of soft robots using model order reduction

Keywords: Soft Robotics, Robot simulation, Model Order Reduction, Proper Orthogonal Decomposition, Hyperreduction, Energy conserving sampling and weighting Method

Obtaining an accurate mechanical model of a soft deformable robot compatible with the computation time imposed by robotic applications is often considered as an unattainable goal. This paper should invert this idea. The proposed methodology offers the possibility to dramatically reduce the size and the online computation time of a Finite Element Model (FEM) of a soft robot. After a set of expensive offline simulations based on the whole model, we apply snapshot-proper orthogonal decomposition to sharply reduce the number of state variables of the soft robot model. To keep the computational efficiency, hyperreduction is used to perform the integration on a reduced domain. The method allows to tune the error during the two main steps of complexity reduction. The method handles external loads (contact, friction, gravity...) with precision as long as they are tested during the offline simulations. The method is validated on two very different examples of FE models of soft robots and on one real soft robot. It enables acceleration factors of more than 100, while saving accuracy, in particular compared to coarsely meshed FE models and provides a generic way to control soft robots.

I. Introduction

Soft Robotics is an emerging field of robotics. Unlike traditional robots typically made of articulated rigid bodies, soft robots are made of soft materials, such as silicone, and their movements are created by deformation. The design of the soft robots is often inspired by nature: elephant's trunk [START_REF] Hannan | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF], octopus' arm [START_REF] Laschi | Soft robot arm inspired by the octopus[END_REF] or worms for example. It is thought that for some applications, soft robots may be more suited than their rigid counterparts. These applications may be manipulation of fragile objects, locomotion in difficult terrain, interaction with humans, or surgical applications to name a few. The potential of soft robots were identified in [START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF], [START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF]. The suitable mechanical impedance provided by the soft materials may also be achieved by Variable Stiffness Actuators (VSAs) or Series Elastic Actuators(SEAs) [START_REF] Della Santina | Controlling soft robots: balancing feedback and feedforward elements[END_REF].

A. Simulation

The simulation of robots is an important field in robotics. Simulators provide assistance for robot prototyp- 1 DEFROST team, Inria Lille-Nord Europe, 40 avenue du Halley, 59000 Lille, France olivier.goury@inria.fr 2 University of Lille -CRIStAL UMR CNRS 9189, 59655 Villeneuve d'Ascq Cedex France ing, testing of control methods, online or offline planning of the motion within the environment, reward evaluation tests for learning methods or supervision. General platforms like V-Rep [START_REF] Freese | Virtual robot experimentation platform v-rep: A versatile 3d robot simulator[END_REF], Gazebo [START_REF] Koenig | Design and use paradigms for gazebo, an open-source multi-robot simulator[END_REF] or even more specialized simulation system, like openHRP [START_REF] Kanehiro | Openhrp: Open architecture humanoid robotics platform[END_REF] dedicated to humanoid robotics, provide mechanical simulation of articulated rigid robot dynamics. Interaction of robots with their environment could also be simulated thanks to collision detection and response available in physics engines like ODE [START_REF] Smith | Open dynamics engine[END_REF], PhysX [START_REF] Rieffel | Evolving soft robotic locomotion in physx[END_REF], Bullet [START_REF] Coumans | Bullet physics engine[END_REF]. See [START_REF] Erez | Simulation tools for modelbased robotics: Comparison of bullet, havok, mujoco, ode and physx[END_REF] for a comparison of these tools. The code of these engines is optimized to obtain very fast simulation, in particular to allow interactive and real-time simulation of the robot, which are very intuitive and useful for the users. However, such physics engines with real-time level of performance are not available for deformable robots in the general case.

In some particular cases, for instance when the robot is considered as a series of deformable rods (like concentric tube robots), very fast and accurate models could be used based on Kirchoff [START_REF] Rucker | geometrically exact model for externally loaded concentric-tube continuum robots[END_REF] or Cosserat rod theory [START_REF] Spillmann | Corde: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects[END_REF], [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF]. To model soft robots, tools such as Voxelyse [START_REF] Hiller | Dynamic Simulation of Soft Multimaterial 3D-Printed Objects[END_REF] were developed, based on a representation of the robot structure in Voxels supported by a lattice modelled using beam theory. For more realistic mechanical simulations, finite element (FE) based method can be used to predict the robots deformation with great accuracy (see for example [START_REF] Zhang | Design and analysis of soft grippers for hand rehabilitation[END_REF]). However that is an expensive method that does not allow a priori for interactivity. Real-time FE simulation and control of soft robots [START_REF] Duriez | Control of Elastic Soft Robots based on Real-Time Finite Element Method[END_REF] were made possible using the SOFA framework 1 , which is dedicated to the fast simulation of deformable bodies, and the SoftRobots plugin 2 . This method was used for various type of robots geometries and actuators [START_REF] Zhang | Kinematic modeling and observer based control of soft robot using real-time finite element method[END_REF], [START_REF] Coevoet | Optimization-based inverse model of soft robots with contact handling[END_REF], [START_REF] Duriez | Framework for online simulation of soft robots with optimization-based inverse model[END_REF]. However, the real-time constraints made the method possible only for relatively coarse meshes and simple material constitutive laws.

The objective of this study is to go much further on this limit of computation time and to allow a real-time simulation of flexible robots on meshes and models much more complex, while remaining generic. For that purpose, we propose to use a technique of model order reduction. This code is developed in a SOFA plugin 3 . In the future, the SOFA framework could be integrated in more traditional simulation frameworks like Gazebo, to provide a physics engine to simulate soft robots.

B. Model order reduction

Model order reduction methods are well-used in the computational mechanics community to reduce the computational cost of FE simulations, with material design or virtual-prototyping applications in mind. Model order reduction (MOR) methods consist in projecting the large FE equations of motion onto attractive subspaces of smaller dimensions allowing much faster computations [START_REF] Sifakis | Fem simulation of 3d deformable solids: A practitioner's guide to theory, discretization and model reduction[END_REF].

In practice it allows to dramatically reduce the number of degrees of freedom of the computed mechanical system, while keeping the precision. One of the main methods is the snapshot-POD (Proper orthogonal decomposition) [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. part I: coherent structures[END_REF], [START_REF] Goury | Computational time savings in multiscale fracture mechanics using model order reduction[END_REF], [START_REF] Carlberg | Efficient non-linear model reduction via a least-squares petrov-galerkin projection and compressive tensor approximations[END_REF], which proposes to represent the highdimensional deformation of the structure within a small attractive subspace of much smaller dimension, defined by a basis of a few vectors. To find this small but representative subspace, the process involves an offline stage where all the potential movements of the structure studied are sampled and stored in what is called the snapshot. This offline stage is computationally intensive, since many fine simulations are performed, but it is performed only once to build the reduced order model. The snapshot space is condensed in a reduced basis using a singular value decomposition and a truncation. The resulting reduced-model of the structure has few degrees of freedom and the simulation can be performed at a much higher rate. To keep the method numerically efficient in the case of nonlinear material behaviour, some hyperreduction method [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF], [START_REF] Carlberg | Efficient non-linear model reduction via a least-squares petrov-galerkin projection and compressive tensor approximations[END_REF], [START_REF] Farhat | Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency[END_REF], [START_REF] Hernandez | Dimensional hyper-reduction of nonlinear finite element models via empirical cubature[END_REF] are used to reduce that computation. Some notable contributions in the idea of using model order reduction technique for fast FE simulations can be found for a simple soft robot example using Proper Generalised Decomposition (PGD) [START_REF] Chenevier | Reduced-order modeling of soft robots[END_REF], or in the context of computer animations [START_REF] Sin | Vega: Non-linear fem deformable object simulator[END_REF], [START_REF] Barbič | Real-time subspace integration for st. venant-kirchhoff deformable models[END_REF] or material design [START_REF] Xu | Interactive material design using model reduction[END_REF]. These latter methods generate the reduced basis by computing vibration modes and their derivatives to account for large deformations.

C. Contribution

The contribution of the paper is, to our best knowledge, the first published method to obtain a very fast mechanical model of a soft robot -compatible with real-time robotic constraints -based on a complex FE model reduced by a snapshot-POD method.

To demonstrate efficiency and genericity, the method will be tested on three very different cases: one deformable cable-driven parallel manipulator, a multigait pneumatic robot for locomotion, showing that the method is compatible with contact constraints with the floor, and a tendondriven soft tentacle interacting with its environment. The method enables very impressive acceleration factors, while keeping a good accuracy and in all cases true real-time (computation time = simulation time) could be achieved. In the third example, the method is used to control a real tentacle soft robot.

Despite some limitations, especially on local deformations that are discussed in the paper, the approach is very comprehensive and opens completely new perspectives in terms of modeling and simulation for deformable robots.

II. Soft robots finite element model

In this section we define the finite element modeling framework used to model the deformation of soft robots in general. The robot may have any shape. Consider the deformation of a soft-robot subjected to external forces and contacts. Its dynamics is defined using Newton's second law:

M(q) v = P(t) -F(q, v) + H T λ, ( 1
)
where q is the vector of position, v is the vector of velocity, M(q) is the mass matrix, P gathers external forces, F accounts for internal forces (depending on the material the soft robot is made of) and H T λ is the vector of constraint forces contributions, which are typically generated when the soft-robot boundary enters in contact with its environment, or when an actuator applies force on the soft-robot.

A. Space discretization

In this paper, we mainly focus on 3D geometries though the method would be applicable to shell models or beam models in principle. The geometry of the robots is meshed with tetrahedrons or hexahedrons. The finite element method is applied to express the deformation of the robot in a discrete way. In this paper, we will use linear tetrahedrons though in principle, the method would work for any order of elements.

B. Implicit time integration

The continuous time space is subdivided into discrete intervals: [t 0 , t 1 , ..., t nt ]. Let's consider a small time interval [t n , t n+1 ], and let h = t n+1 -t n . Integrating equation (1) over this interval gives (assuming the mass matrix is constant):

M(v tn+1 -v tn ) = tn+1 tn (P(t) -F(q, v)) dt + hH T λ (2) q tn+1 -q tn = tn+1 tn v(t)dt (3)
Using an implicit Euler integration scheme, we obtain:

M(v tn+1 -v tn ) = h p tn+1 -F(q tn+1 , v tn+1 ) + hH T λ (4) q tn+1 = q tn + hv tn+1 , (5)
where p tn+1 is the value of the function P at time t tn+1 (known a priori). The internal forces F is a nonlinear function of the positions and the velocities, which means equation ( 4) can not be solved exactly in one step and one can use an iterative method such as Newton-Raphson.

A less expensive method, though less accurate when the timestep is large compared to the internal forces nonlinearity, is to solve the system by applying a Taylor series expansion to F leading to the following first order approximation:

F(q tn+1 , v tn+1 )) = F (q tn + dq, v tn + dv) (6) = f tn + δF δq dq + δF δv dv (7)
Using dq = q tn+1 -q tn = hv tn+1 and dv = v tn+1 -v tn , we obtain:

M + h δF δv + h 2 δF δq A(qt n ,vt n ) dv (8) =-h 2 δF δq v tn -h f tn + p tn+1 b(qt n ,vt n ) +hH T λ. (9)
This is exactly equivalent to perform the first iteration of the Newton-Raphson algorithm. In this paper, we choose this linearisation strategy for performance purposes, but the method we will describe in the following section would also work using the Newton-Raphson method. The computational cost of one Newton-Raphson iteration is identical to one linearisation step. When the mesh used to describe the soft robot is fine, equation ( 8) is of large dimension and prohibits an interactive simulation of the soft robot deformations.

C. Constitutive law of the soft robot material

The constitutive law of silicone is generally recognised as being hyperelastic, which establishes a nonlinear relationship between stresses and strains in the material. In this paper, to model the behaviour of the soft robot, we will consider two kinds of constitutive laws:

• Co-rotational linear elastic. This formulation assumes a linear stress-strain relationship within the frame attached to each element of the FE mesh. The corotational element stiffness contribution is computed using the rotation transformation between the deformed element compared to its rest position. This results in a non-linear constitutive law, since the displacements are needed to compute the stiffness matrix. This law can actually represent large deformations, as long as those deformations involve only large rotations but relatively small deformations of each individual elements. • Mooney-Rivlin Hyperelastic. This is a well known model for hyperelasticity [START_REF] Mooney | A theory of large elastic deformation[END_REF]. It assumes the following strain energy density function:

W = C 01 ( Ī1 -3) + C 10 ( Ī2 -3) + D(J -1) 2 , ( 10
)
with Ī1 = J -2/3 I 1 , Ī2 = J -4/3 I 2 , J = det(F), F being the deformation gradient, I 1 and I 2 being the first and second strain invariants of the right Cauchy-Green deformation tensor.

These two constitutive laws allow to compute the internal forces vector F through the integration of partial differential equations coming from continuum mechanics.

To have a good accuracy on the numerical model, it is often necessary to use a very fine discretisation which leads to systems of high dimension. Hence, there is a real need for reduction.

III. Reduction of the full order Finite element model

In this section, we present a general framework for the reduction of the full order dynamic model describing the robot's motion presented in the previous section. This reduction will be performed in two stages: a projection-based method to express the state variables in a low-dimensional subspace, and a system approximation method, to integrate the mechanical properties only on a small subset of elements of the large FE mesh.

A. Snapshot -Proper Orthogonal Decomposition

The proper orthogonal decomposition (POD) [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. part I: coherent structures[END_REF], [START_REF] Carlberg | Efficient non-linear model reduction via a least-squares petrov-galerkin projection and compressive tensor approximations[END_REF], [START_REF] Goury | Computational time savings in multiscale fracture mechanics using model order reduction[END_REF] is a method used widely in the computational mechanics community to generate relevant bases to be used in reduced-order models. The model to reduce is typically parametrized, and the reduced model is then valid for any values of the parameter. To apply the method in our context, we define Λ, the parameter space of scenarios of constraints applied onto the robot (that may be from actuators or external contacts) parametrized by the constraints 4 λ.

Then, the POD decomposes the position q(t), solution of (8) in a truncated expansion of orthonormal vectors (where we show clearly the dependency of q on the constraints λ):

q (t, λ) ≈ q(0) + N i=1 φ i α i (t, λ) (11)
= q(0)

+ Φ α(t, λ), (12)
such that the orthonormal basis (φ 1 , φ 2 , . . . , φ N ) is minimizing the cost function representing the sum of all the projection errors of the exact solution onto the basis:

J(Φ) 2 = λ * ∈Λ t=tn t t=t0 q(t, λ * ) - φ i T q(t, λ * ) φ i 2 2 dt. (13)
In practice, this technique is not useful since it requires to know the exact solution at each point of time and for every scenario of loading. The snapshot-POD proposes to use a discrete approximation of this cost function, where we define a discrete subset of Λ: Λ = {λ 0 , λ 1 , . . . , λ n λ }. The cost function becomes:

J(Φ) 2 = λ * ∈ Λ t=tn t t=t0 q(t, λ * ) - φ i T q(t, λ * ) φ i 2 2 . ( 14
)
The set q(t, λ)| t ∈ [t 0 , t nt ], λ ∈ Λ can be assembled in a snapshot matrix S = [q(t 0 , λ 0 ), . . . , q(t nt , λ n λ )]. The solution to the problem of minimising the cost function defined in [START_REF] Spillmann | Corde: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects[END_REF], under the constraint of orthonormality of the basis is given by performing a singular value decomposition (SVD): S = UΣV T , with U and V unitary matrices and Σ a rectangular diagonal matrix containing the singular values of S. The optimal basis of order p is then given by the first p left singular vectors (in U) associated to the p largest singular values σ i . The truncation error of a POD transform of order p is given by:

ν 2 = nq i=p+1 σ 2 i nq i=1 σ 2 i (15)
The reduced basis Φ, minimising the cost function ( 14), will give a good approximation of the position q if the set Λ is exhaustive enough and represents well all the workspace of the robot (induced by all possible actuations) and the contacts the robot may encounter. Note that computing the cost function J requires the knowledge of the full order solution over the set Λ, and it is a very expensive step. However, it is performed only once, and the reduced basis Φ is computed once for all. The definition of the set Λ, which is the set of scenarios of loading the robot will be subjected to in the offline phase, has a strong influence on what deformations the resulting reduced order model will then be able to reproduce. In the context of this paper, this set is defined by all possible combination of actuation. Assuming there are N A actuators, and defining a range divided in a certain number of samples s i for each actuator, the number of possible actuations is Π N A i=1 s i . For example, a robot with 4 actuators with each actuator having a range divided in 10 steps leads to 10 4 = 10000 configurations. In this paper, we will simplify the generation of the snapshot by only considering the 2 extreme values of the range of each actuator, leading to 2 N A configurations and the robot will be simulated while moving directly from one actuation state to another, as described in Figure 1. Despite this simplification, the number of tests to perform to generate the snapshot grows exponentially with the number of actuators. In the case the number is not tractable with the computational power at hand, non-exhaustive sets can be selected through Latin-Hypercube sampling. In the general context of model order reduction, the problem of selecting an optimal snapshot to create accurate reduced models is an active research field. Beyond uniform sampling, other strategies exist to sample the workspace: random selection [START_REF] Goury | Automatised selection of load paths to construct reduced-order models in computational damage micromechanics: from dissipation-driven random selection to bayesian optimization[END_REF], error estimation to select the samples in an iterative way [START_REF] Quarteroni | Certified reduced basis approximation for parametrized partial differential equations and applications[END_REF], gradient-optimisation [START_REF] Bui-Thanh | Model reduction for large-scale systems with high-dimensional parametric input space[END_REF] or statistical inference [START_REF] Paul-Dubois-Taine | An adaptive and efficient greedy procedure for the optimal training of parametric reduced-order models[END_REF], [START_REF] Goury | Automatised selection of load paths to construct reduced-order models in computational damage micromechanics: from dissipation-driven random selection to bayesian optimization[END_REF]. These strategies may be used in future works. Fig. 1: Example of actuation sequence for a robot with 2 actuators and hence 2 2 = 4 different extreme configurations. In 3 steps, the actuation moves from one configuration to another. This sampling of the space of actuation is not exhaustive but tends to give sufficient information to build an accurate reduced order model as it explores the extremities of the actuation space.

B. Galerkin projection of the equations of motion onto the reduced space

Once the basis is defined, the reduced expression of the state variables (position, velocity) can be introduced into the equation of motion [START_REF] Kanehiro | Openhrp: Open architecture humanoid robotics platform[END_REF]. Note that by differentiating with respect to time, the reduced expansion [START_REF] Coumans | Bullet physics engine[END_REF] naturally leads to a reduced expression for the vector of velocities: v(t) = q(t) = Φ α(t). Plugging these expressions into equation ( 8) and projecting onto the reduced subspace (thus a Galerkin projection) leads to:

Φ T A(q tn , v tn )Φ Ar ḋ α(t) = Φ T b(q tn , v tn ) br + (H Φ) T Hr λ, (16)
where dα = α tn+1 -α tn and ḋ α = αtn+1 -αtn . Note that we performed a Galerkin projection (using the reduced basis itself as projector) for simplicity, but a Petrov-Galerkin projection could have been used (a study in the two different approaches was made in [START_REF] Carlberg | Efficient non-linear model reduction via a least-squares petrov-galerkin projection and compressive tensor approximations[END_REF]). This is a reduced equation that is of the dimension N of the subspace generated by the basis Φ. The new unknown becomes dα which is of much smaller dimension than the vector of positions q and of velocities v.

C. System approximation of the Finite Element mesh for fast efficient construction of the nonlinear operators

The reduced equation ( 16), though of smaller dimension than the original full order model, is still expensive to build. Indeed, because the equations of motion are nonlinear, the operator A r is dependent on the current state of the system and can not be pre-computed. Hence, it is still necessary to loop over all the elements of the high dimensional mesh to evaluate the full order operator A before projecting it onto the reduced basis Φ. The same is true for the right hand side b r . As a result, the gain in computational time provided by the reduced model will be negligible. This issue is well-known in the field of reduced-order modeling and can be tackled by applying a second approximation called hyperreduction or system approximation. This method provides a way to solve the reduced equation while computing the mechanical properties on a small subset of the elements of the full order mesh only. The procedure has several variations in the literature. We can cite the hyperreduction method [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF], the discrete empirical interpolation method [START_REF] Chaturantabut | Discrete empirical interpolation for nonlinear model reduction[END_REF] or more recently the energy-conserving sampling and weighting (ECSW) method [START_REF] Farhat | Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency[END_REF] or empirical cubature method [START_REF] Hernandez | Dimensional hyper-reduction of nonlinear finite element models via empirical cubature[END_REF]. In this work, we focus on the ECSW method (very similar to the empirical cubature method) which builds stable and consistent reduced models. First, for the unique sake of simplicity of the notations, we rewrite equation ( 16):

Φ T g(t n ; λ) = (H Φ) T λ, (17)
with g(t n ; λ) = A(q tn , v tn )Φ ḋ α(t) -b(q tn , v tn ). Here, we show clearly the dependence of the operator g on the constraints λ, which can be seen as a parameter of the simulation. The method proposes to go back to the element level. Indeed, the left-hand side can be decomposed as:

Φ T g(t n ; λ) = Φ T e∈E Y (e) g e (t n ; λ) (18) = e∈E Φ T e g e (t n ; λ), (19)
where Y (e) ∈ R nq×ne is the assembly operator that maps degrees of freedom from the element level (e) to the global degrees of freedom of the entire system. Y (e) is rectangular and only filled with zeros and ones. Hence, Φ e is the restriction of Φ to the degrees of freedom corresponding to element (e). g e (t n ; λ) is the contribution of element (e) to the global term g(t n ; λ). Then, the ECSW method interprets the quantities Φ T e g e (t n ; λ) as a set of energies and therefore proposes to approximate g by performing the sum over only a small subset of the elements E (the reduced integration domain or RID) weighted with some positive coefficients ξ e :

Φ T g(t n ; λ) ≈ e∈ E ξ e Φ T e g e (t n ; λ). (20)
This new expression is much less computationally intensive as the computation of all the operators will involve a loop over a small number of elements rather than the entire high-dimensional mesh. Now, the set of elements of the reduced integration domain E and the corresponding weights (ξ e ) e∈ E can be found through an optimization problem, based on a database of unassembled elements contributions computed previously (This implies that another set of offline simulations needs to be performed to store these element contributions, and typically, these are done using the same actuation sequence that generated the snapshot in the previous stage). The problem can be expressed in the following way:

ξ = arg min ξ * 0 ( 21
)
subject to:

Gξ * -y 2 ≤ τ y 2 ξ * ≥ 0 ,
where τ is a tolerance defined by the user. The zero norm x 0 is the number of non-zeros entries in x. G = [G ke ] k∈ Λ,e∈ E is a block matrix containing the unassembled element contributions for each configuration λ ∈ Λ, with the column block G ke = Φ T e g e (t k , λ k ) for the k th timeparameter configuration and element number e. y = [y k ] k∈ Λ is the assembled counterpart:

y k = Φ T g(t k , λ k ) = e∈ E G ke .
This minimisation problem is a non-negative least square problem [START_REF] Lawson | Solving least squares problems[END_REF] and finding the optimal solution typically involves an intractable exhaustive combinatoric search. In practice, a suboptimal greedy algorithm, like described in [START_REF] Peharz | Sparse nonnegative matrix factorization with 0 -constraints[END_REF], [START_REF] Farhat | Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency[END_REF]), can be used to find an acceptable solution to problem [START_REF] Duriez | Framework for online simulation of soft robots with optimization-based inverse model[END_REF].

D. Summary for building the reduced order model of a soft robot

In this section, we sum up the main steps of the reduction workflow in the form of an algorithm: Algorithm 1 Workflow for building reduced order model of a soft robot Ensure: The accuracy of the full order model is representative of the real soft robot.

-Evaluate number of actuators N A and their range of actuation.

-Define number of samples s i within each actuator range (minimum is 2) if it is feasible to compute Π N A i=1 s i simulations to generate an exhaustive snapshot then Compute and store the positions (or state variables in general) in the snapshot else Define an alternative non-exhaustive snapshot selection strategy (Latin Hypercube, Bayesian, etc...) to compute the snapshot end if -Perform a singular value decomposition on the snapshot matrix and truncate up to a tolerance ν defined in [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF]. This defines the reduced basis Φ.

-Compute the reduced model on the snapshot space again, storing the FE elements contributions along the way.

-Define a tolerance τ and solve the associated NNLS problem [START_REF] Duriez | Framework for online simulation of soft robots with optimization-based inverse model[END_REF]. This gives the reduced integration domain and the weights connected to it.

IV. Practical applications

The proposed method allows to reduce greatly the computational time of an FE model of a soft robot. First, two examples on two very different kinds of robots will be presented to show the method can be generic. The practical applications are hence primarily concentrated on the numerical simulation, assuming that a converged FE solution on a fine mesh can be used as "ground truth". In those sections, we show the application of the method on these two robots and analyse performance in terms of accuracy and computational time with respect to the full order model. Then, the applicability of the method will be shown to control a real tentacle robot.

A. A first example: fixed soft silicone robot actuated by four cables

Let us consider the parallel soft silicone robot described in [START_REF] Duriez | Control of Elastic Soft Robots based on Real-Time Finite Element Method[END_REF]. This purely academic robot is displayed in Figure 2. It is entirely made of soft silicone and is actuated by four cables controlled by step motors located at its center. Pulling on the cables has the effect of lifting the effector located on top of the robot. The "game" with this robot is to control the position of the effector by pulling on the cables. Fig. 2: From left to right: coarse mesh with 1400 nodes; fine mesh with 15553 nodes; reduced mesh used for hyperreduction : 861 nodes and only 376 tetrahedra. a) Finite element mesh refinement: In [START_REF] Duriez | Control of Elastic Soft Robots based on Real-Time Finite Element Method[END_REF], the robot was controlled through real-time finite element simulation based on a mesh of 1628 nodes and 4147 tetrahedra. That size of mesh was manageable in real-time on a standard desktop computer. The simulation made using this underlying mesh was accurate enough to control the robot, only considering the displacement of the effector point, located on the top of the robot and with a limited range on the pulling of the cable actuators. However, even the full convergence of the FE model could require a larger number of elements (convergence of the displacement or of the stress field inside the domain which could be important for the equilibrium with external loads). In some configurations, in particular with large deformations, the actual position each of the four arms of the robot may not be accurately predicted for example. When considering an application where the robot arms may enter in contact with the environment, an accurate prediction of their position becomes relevant. In Figure 3, we show the difference in prediction between a coarse mesh and a much finer one. In this case, the error in displacement is concentrated in the region of the "elbows" of the robot. To get an idea of how fine the mesh should be to be accurate enough, we first perform a quick sensitivity analysis on the behaviour of the robot with increasingly finer meshes. We compare the results using four different meshes of increasing size. All the meshes are generated with the CGAL mesh generator 5 . Results are shown in Figure 3. From these experiments, we decide to pick the mesh with 15553 nodes, which presents a relative error of less than one percent with the finest mesh. Indeed, typically, the reduced order model will generate an error of that magnitude. It is hence unnecessary to go further in mesh refinement, unless a much greater accuracy is needed. Fig. 3: (Left) Superposition of the robot deformation using two meshes with different refinement in one scenario of cables pulling. We can see that the differences between the fine model and the coarse model of the robot mostly concentrates along the arms. (Right) Relative error in the displacement of the arms of the robot with increasingly finer meshes. The error is computed relative to a finer FE mesh with 34118 nodes and 171603 elements. We consider that the 3 rd mesh with 15553 nodes is fine enough, given that the accuracy barely improves with a mesh more than twice as fine.

b) Finite Element model used: We will consider two full FE models using two different constitutive laws:

• A linear elastic law with a co-rotational formulation. The young Modulus is set to 450 MPa, which corresponds to sensible value for silicone. • A hyperelastic law, namely a Mooney-Rivlin law. For consistency, the parameters values C 01 , C 01 and D (from ( 10)) are chosen to match those of the linear elastic law in the limit of small displacements. c) Generation of the reduced basis: We proceed by "shaking" the robot following each cable actuator range. A range of [0, 40]mm is defined for each actuator, and the snapshot phase test all possible combinations of actuating cases at the ends of the range. With these 4 actuators with a 1-dimensional range, this results in a manageable set of simulations to perform. Indeed, we want to make the robot reach all 2 4 = 16 different extremity states that are accessible within the range of the actuators. Dividing the loading from one state to another in 9 steps, this provides a set of 9 * 15 + 1 = 136 positions of the robot, following the same principle displayed in Figure 1.

Performing a singular value decomposition to obtain the reduced basis, we define an error truncation criteria ν = 10 -3 . From equation [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF] this criteria can be achieved with a selection of 34 basis vectors when using the elastic law, and 30 when using the Mooney-Rivlin hyperelastic formulation. This error criterion is only valid on the snapshot space. If that snapshot space is exhaustive enough, the actual error between the true model and the reduced model should be quite small. Some of those vectors are displayed in Figure 4.

The system approximation method (ECSW), described in the previous section, is then applied to reduce the computational costs with a tolerance τ varying from 0.03 to 0.1 in the minimisation [START_REF] Duriez | Framework for online simulation of soft robots with optimization-based inverse model[END_REF]. This leads to reduced integration domains of various sizes.

d) Numerical Results: In Tables I andII, we compare the accuracy of the simulation using the coarse mesh versus using the reduced order model with various sizes. A more graphical representation of the performances of the reduced order model is displayed in Figure 6. The error is measured with respect to the full order model evaluated on the fine mesh. Several remarks can be made.

• First of all, we can see that the reduced model can be evaluated at a much faster rate than its full order counterpart, no matter which constitutive law is used (elastic or hyperelastic). Indeed, for the ROM with 34 basis vectors and a reduced integration domain of 376 elements, the speedup is around 58. • As expected, the accuracy is directly linked to the size of the ROM. Fewer elements in the ROM leads to a larger error. • It becomes clear, especially on Figure 6, that, with an appropriate size of the ROM, a frame rate at least as large as the full order model on the coarse mesh model can be achieved while having a much better accuracy. In other words, the ROM based on the expensive but accurate full order model on the fine mesh keeps a good accuracy while being computationally as inexpensive as a full order model based on a much coarser mesh and hence quite inaccurate. This is particularly true for the simulation using the hyperelastic formulation, where the accuracy is almost an order of magnitude better than the full coarse model while having a runtime twice as low. In the context of simulation of soft robots for realtime control applications, this opens up the possibility to simulate interactively the behaviour of soft robots with a higher accuracy. In particular, this allows to start from a real accurate digital description of the robot with the finite element method before applying the reduced order modelling workflow to make the simulation computationally tractable within the real-time control applications constraints.

B. More complex robot geometry and predictable contact handling: Multigait Soft Robot

In this example, we consider the multigait soft robot, as described in [START_REF] Shepherd | Multigait soft robot[END_REF]. The robot is displayed in Figure 7. This robot is made of two layers: one thick layer of soft silicone Relative errors between the simulation with the coarse mesh and the reduced order models and the finer mesh simulation using the Mooney-Rivlin Hyperelastic constitutive law. In that case, the hyperreduced model leads to a much smaller error than the full order model with the coarse mesh while still running about 2 times faster.

containing the cavities, and one stiffer and thiner layer of Polydimethylsiloxane (PDMS) that can bend easily but does not elongate. The robot is actuated by five air cavities that can be actuated independently. The effect of inflating each cavity is to create a motion of bending. Then, by actuating with various sequences each cavities, the robot can move along the floor. We propose here to make an interactive simulation of the robot crawling on the floor. a) Finite Element model of the multigait soft robot: To model the multigait soft robot, we use two different internal force contributions, one volumetric one for the main body of the robot, and one surfacic one for the nonstretchy layer of PDMS. A membrane model is used for that layer: bending forces are ignored, and internal forces are only created by stretch. Using this 2D representation, we can avoid meshing that thin layer in 3D which would have required many elements. This decomposition is displayed in Figure 8. For the main body, internal forces are computed as a 3-dimensional force field, ans a corotational elastic constitutive law is used. Young's modulus constitutive parameter is taken from literature, to match the material used in [START_REF] Shepherd | Multigait soft robot[END_REF]. For the material properties of the membrane representing the PDMS layer, a much higher stiffness is used to model the fact it is not extensible.

To model the actuation due to the pneumatic cavities, we do not use a dynamic model for the pressurised air and its interaction with the cavities, but rather assume a uniform pressure dispatched regularly over the cavities in a direction orthogonal to the surface of the cavity. This approach is explained in more details in [START_REF] Duriez | Framework for online simulation of soft robots with optimization-based inverse model[END_REF]. In Figure 9, we can see the effect of the PDMS layer on the deformation with several meshes and reduced order models (ROM) with increasing reduced integration domains, in terms of accuracy versus computational time. The error is measured on the goal point and on the arms of the robot (for the hyperelastic model, we display the error for the arms only for clarity of the graph). Values are extracted from Tables I and II. For the FOMs (continous lines), a more refined mesh implies a better accuracy, but at a high computational cost which results in a timestep rate decrease. To reach Timestep rates allowing interactivity (circa 25 fps for example), the FOM has to be based on a coarse mesh and hence has low accuracy. Comparatively, for the ROMs (discontinuous lines), the interactivity rate can be achieved with an error of an order of magnitude smaller than the FOMs. The timestep rate can increase dramatically while decreasing the size of the RID but the error remains low. Fig. 7: Multigait soft robot. Picture taken from [START_REF] Shepherd | Multigait soft robot[END_REF] of the legs while inflating the corresponding cavity. The layer has the clear effect of bending the leg, which matches the behaviour of the real robot. In Figure 10, we can see a comparison between the real robot and its finite element simulation.

Fig. 8: Zoom in the mesh: a membrane constitutive law is applied onto the bottom layer of the soft robot. That contribution is added to the nodes located at the bottom of the robot. b) Generation of the reduced basis and rigid body modes: Unlike the previous example, this robot is not fixed on the floor (so a dynamical model is necessary for the simulation) , and the contacts of the floor onto the robot influences its deformation. To generate a relevant snapshot, we perform the offline phase while simulating the robot inside its environment, i.e. onto the floor, to account for those specific deformations but without friction. So during the offline phase, the robot is not moving, it is only deforming while staying at its initial position. The robot may also crawl around on the floor since it is not Fig. 9: Effect of the bottom layer on the deformation of the robot captured by the simulation. Left: without bottom layer, the material expands in no special direction. Right: with the bottom layer, the material does bend. Fig. 10: Comparison between a fine FE model of the soft robot and the real robot (from [START_REF] Shepherd | Multigait soft robot[END_REF]) with identical pressures in the cavities. The finer mesh is made of 15787 nodes and 59296 elements. In that case, we can see that it reproduces the behaviour of the real robot fairly well. fixed. To allow the robot to move, the reduced basis should include the rigid body motions R (i.e. translations and rotations) of the robot to be able to move freely onto the floor. The reduced basis computed from the snapshot does not include those rigid modes R, corresponding to pure translations (R t ) and pure rotations (R θ ), which do not create internal forces on the robot. Hence, it is necessary to concatenate them to the reduced basis Φ, to allow for these movements on the reduced model, while erasing their projection onto it to keep the orthogonality of the basis and a good conditioning of the reduced system:

Φ = [R , Φ snap - i∈{Rt,R θ } R T i Φ snap R i ], (22)
where Φ snap is the reduced basis computed with the POD from the snapshot set. The robot has 5 pneumatic actuators, which means the snapshot generation will involve at least 2 5 = 32 extreme configurations, if considering only two cases for each actuator (no air pressure and maximum air pressure). For the sake of stability of the simulation, we do not apply the pressure all at once but iteratively between each configuration, over 15 steps. This leads to a snapshot with about 500 configurations. Performing the singular value decomposition with an error tolerance ν = 10 -3 (respectively 10 -2 ) from equation [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF] leads to a basis of 60 modes (respectively 30). Adding the rigid body modes for translation (we neglect the rotations) leads to a basis with 63 vectors. In Figure 11 efficiency, we proceed to apply a method of hyperreduction. The application is slightly different in this case, since there are two distinct internal forces applying on the robot, a three-dimensional one on the body, and a two-dimensional one on the bottom surface. The ECSW procedure is applied independently for each internal force to output a reduced integration domain for each force. The resulting domains are displayed in Figure 12. d) Speed-up and behaviour of the reduced-order model: Computational time gains are displayed in Table III. We can see that the reduced model is up to more than 100 times faster than the full order model. One interesting point is that unlike the previous example of the parallel cable robot, the reduced model without hyperreduction already has some noticeable gains in term of computational time over the full order model. This is due to the fact that contact handling can be computed in a much faster way with the reduced model since the sytems involved are of much smaller dimensions. In the video material associated with this paper, it can be seen that the reduced simulation of the robot is able to reproduce the undulating movement of the robot. In particular, it is able to replicate the forward movement of the robot following the same sequence of actuation.

C. Application to the control of a real cable-driven soft robot using an inverse model with contact handling

To show the applicability of the method to inverse kinematics and control of a real robot, we apply the procedure to the tentacle robot described in [START_REF] Coevoet | Optimization-based inverse model of soft robots with contact handling[END_REF]. This robot is divided in two sections, and actuated by 8 cables, 4 for each section. It also has one actuator allowing to move the whole robot forward in a translation. A drawing of its geometry is displayed in Figure 14. In this section, we will control the robot through inverse kinematics. A endeffector is defined at the tip section of the tentacle robot, and an inverse model based on the reduced order finite element simulation computes the actuation needed for the end-effector to reach a user-defined goal. The inverse problem formulation minimises the distance between the goal and the end-effector, while having a penalty for the mechanical energy of the robot deformation, which guarantees the uniqueness of the solution at each time step. The optimisation problem is solved using a Quadratic program with linear complementarity constraint (QPCC). More details on the definition of the inverse problem as well as its solving can be found in [START_REF] Coevoet | Optimization-based inverse model of soft robots with contact handling[END_REF]. We define a fine discretisation with 28221 tetrahedra (compared to a coarser description with 1565 tetrahedra in [START_REF] Coevoet | Optimization-based inverse model of soft robots with contact handling[END_REF]). Applying the the reduced order modelling strategies, the robot is tested in 2 8 = 256 extremity cases of its workspace. With a tolerance set to 10 -3 , this leads to a reduced model of dimension 33, and with a tolerance τ = 0.07, to a reduced integration domain with 102 elements. The meshes are displayed in Figure 15. Though the coarse mesh was sufficient to control the tentacle robot in [START_REF] Coevoet | Optimization-based inverse model of soft robots with contact handling[END_REF] with inverse kinematics, the reduced model allows for a greater accuracy, as can be seen in Figure 16. The reduced model simulation is then tested to control the real tentacle robot using an inverse model with contact handling, in the same way this was done in [START_REF] Coevoet | Optimization-based inverse model of soft robots with contact handling[END_REF]. See Figure 17. The reduced model is compatible with the inverse method for control already used in [START_REF] Duriez | Control of Elastic Soft Robots based on Real-Time Finite Element Method[END_REF], [START_REF] Coevoet | Optimization-based inverse model of soft robots with contact handling[END_REF], and can also handle contacts in this case. Indeed, the deformation created by the contact of the tube onto the tentacle robot remains within the deformation space explored during the offline phase to build the reduced model.

V. Discussion

Applying model order reduction allows to have an accurate mechanical model of the robot, typically of the same order than an FE model of 50000 elements, and to achieve a real-time simulation constraint. It makes possible the use of the dynamical model of the robot for simulation as well as for real-time control (in the parallel robot example and on the tentacle robot, we are testing the inverse model, in the multigait robot example, the dynamic model is tested). The method can be used for contacts with the environment, as long as those contacts are taken into account in the snapshot computation phase (the phase where we shake the robot), or if those contacts do not create deformations outside of that space, like in the tentacle robot example. However, it is not able to cope with local contacts, especially if those where not simulated in the "shaking" phase and create deformations that are unreachable using the actuators only. Another point is that it requires a heavy offline stage which grows exponentially with the number of actuators. The use of our method in the design phase allows to test a virtual representation of the robot at fast rates (with an accurate mechanical model). It allows for instance to design the motion sequence of the actuators or to evaluate the obtained deformations, configurations, workspace. We could also use it to design the control loop. However, at each new change of the design (shape of the robot or type of material used), the offline stage needs to be recomputed. Even if this offline stage takes some hours, it takes much less time (and less efforts) than building the actual robot. In future works, we could avoid some re-computation. For example, beyond the range of actuation, it could be possible to vary the material properties to include them in an interactive design phase. However, this implies that the testing phase becomes even larger. Design involving change in the robot geometry can also be included, especially if the robot is to be made with fundamental building bricks. Reduced models of those buildings bricks can be constructed offline, and added interactively in the design phase. For more general arbitrary geometry changes, the method would require to discretise those geometry changes by creating several reduced models, one for each geometry, but at a tremendous computational cost and is left to future works to optimise.

VI. Conclusion and perspectives

The reduced order modelling technique for soft-robots presented in this paper is very useful to have a mechanical model in real-time and to be able to manage the Fig. 13: Comparison between the hyperreduced model of the soft robot and a full order on a coarser mesh. The coarse mesh is made of 2509 nodes and 8222 elements. We can see that the deformation does not have a great amplitude. This is a known fact that FE methods applied on unrealistically coarse meshes tend to artificially stiffen the structure. On the other hand, the ROM has a good behaviour. In that case, we can see that it reproduces the behaviour of the real robot fairly well. Fig. 15: Finite element meshes of the tentacle: the coarse mesh was used in [START_REF] Coevoet | Optimization-based inverse model of soft robots with contact handling[END_REF] with 1565 tetrahedra, the fine mesh used to build the reduced model with 28221 tetrahedra, the reduced integration domain used for the real-time control with only 102 tetrahedra. compromise between accuracy and computational time. This can be thought of as some sort of machine learning technique, in the sense that a lot of tests are made initially to accumulate data about the robot which is then used to build a reduced model. However, one key difference is that in the case of the reduced-order model, a mechanical model remains, which allows to keep coherency with respect to the physics. This is shown in the case of the multigait soft robot which moves forward while undulating in the reduced-order simulation even if it was tested without friction, and hence remained in the same position during the testing phase. It is however able to move forward in the online phase. Similarly, the snapshot space for the tentacle robot was generated without considering contacts, but the tentacle can deal with contacts in the online stage anyway. Fig. 16: Superposition of the simulation of the tentacle robot using the full order model (FOM) on the coarse mesh, the fine mesh and the reduced order model (ROM). For all three models, the actuation is managed so that endeffector (corresponding to the tip section of the tentacle) reaches a goal materialised by a grey sphere. We can see that the coarse mesh simulation, being artificially stiffer has some error compared to the fine one. The ROM approximates the fine model with a small error.

COMOROS project from which this work is a part of.

Position of the goal controlled by the user

Effector on the real robot

Fig. 4 :Fig. 5 :

 45 Fig. 4: First few basis vectors generated with the snapshot-POD method.

  Linear elastic FOM Goal errors Linear elastic FOM Arms errors Linear elastic ROM Goal errors Linear elastic ROM Arms errors Hyperelastic FOM Arms Error Hyperelastic FOM Arms Error FEM Mesh size increases RID size increases

Fig. 6 :

 6 Fig.6: Comparison between full order models (FOM) with several meshes and reduced order models (ROM) with increasing reduced integration domains, in terms of accuracy versus computational time. The error is measured on the goal point and on the arms of the robot (for the hyperelastic model, we display the error for the arms only for clarity of the graph). Values are extracted from Tables I and II. For the FOMs (continous lines), a more refined mesh implies a better accuracy, but at a high computational cost which results in a timestep rate decrease. To reach Timestep rates allowing interactivity (circa 25 fps for example), the FOM has to be based on a coarse mesh and hence has low accuracy. Comparatively, for the ROMs (discontinuous lines), the interactivity rate can be achieved with an error of an order of magnitude smaller than the FOMs. The timestep rate can increase dramatically while decreasing the size of the RID but the error remains low.

Fig. 11 :

 11 Fig. 11: First basis vectors of the quadruped reduced model. This does not include the translation vectors.

Fig. 12 :

 12 Fig.12: Reduced integration domains for the multigait robot. An independent domain is selected for the three-dimensional force field on one hand, and the twodimensional membrane force field on the other hand. The combination of both leads to the global reduced integration domain for the reduced model of the whole robot.

Fig. 14 :

 14 Fig. 14: Slice view and side view of the tentacle robot.

Fig. 17 :

 17 Fig. 17: Real-time control of the tentacle robot using the reduced order simulation and the inverse model with contact handling.

TABLE I :

 I Relative errors between the simulation with the full order models (FOM) on both coarse and fine meshes and the reduced order models (ROM(size of reduced basis)) and hyperreduced model (HROM(size of reduced basis, size of reduced integration domain)) using a corotational elastic constitutive law.

	Model	Goal Error Arms Error	Tps
	Coarse FOM	0.036	0.048	22
	HROM (30 , 227)	0.0060	0.0049	42
	ROM (30)	0.0050	0.0044	1.0
	Fine FOM	0.000	0.000	0.6

TABLE II :

 II

  are displayed only a few of them.c) Application of the hyperreduction technique: In the same way to the previous example, to get numerical

	Model	Timesteps per second
	HROM 33	20
	HROM 63	12
	ROM 63	0.6
	Fine FOM	0.2

TABLE III :

 III Timestep rate for the full order simulation of the multigait soft robot, the reduced model (ROM), and the hyperreduced model (HROM) with 33 modes and 63 modes.

www.sofa-framework.org

https://project.inria.fr/softrobot

https://project.inria.fr/modelorderreduction

The constraints are actually also defined by the direction parametrized by H but for the sake of simplicity of the notations we only show the constraint dependency with λ instead of H T λ.

www.cgal.org

Acknowledgment

The authors would like to thank the FEDER (European Regional Development Fund ) that co-funded the


Sandra Maria

Neiva Sampaio

Alfredo Kingo

Anne-Elisabeth Elisabeth Laques

Izildinha De

Souza Miranda

Danielle Mitja

Sandra Maria Neiva Sampaio

Oyama Homma

Izildinha Souza

Miranda Souza Miranda

FORMATION OF BABASSU IN DEGRADED AREAS BY HUMAN ACTION IN SOUTHEAST PARÁ, BRAZIL: ITS IMPORTANCE FOR RECOVERY OF LEGAL RESERVE AREAS

Keywords: Babassu palm, Babassu Forest, Legal Reserve Areas, Landscape types, Pastures

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Currently the environmental question in areas of the pioneer front of the Southeastern Pará, Eastern Amazonia has become a great challenge in the search for viable alternatives to conciliate natural resource protection and conservation with sustainable development.

Since the 1960s the occupation process in this region has been fomented by development policies through large projects. In the 1970s the National Integration Program (NIP) influenced the opening of highways and large hydroelectric projects and mineral exploitation [START_REF] Homma | Cronologia da ocupação e destruição dos castanhais no sudeste paraense[END_REF]. This infrastructure supported new migration flows (Homma 2003) and the occupation process culminated with deforestation of large areas for implantation of livestock systems as the basis of economic sustainment of families [START_REF] Pasquis | Rforme agraire" en Amazonie: bilan et perspectives[END_REF][START_REF] Piketty | Les déterminants de l'expansion de l'levage bovain en Amazonie Orientale: consquences pour les politiques publiques[END_REF].

In the case of Itupiranga-PA municipality, Incra data total of settlements corresponds to 5148 families in a total area of 368,197 ha, equivalent to 46.72% of municipal area, where Benfica Settlement Project (Benfica SP) is included. In this settlement, landscape do not represent a homogeneous whole, but involves differentiations that portray in an integrated way, a mosaic of anthropological and natural features, referent to the forest remnants, small areas of family agriculture, large pasture areas and babassu palm forests. This simplification of landscape has made original plant structure unstable, from the point of view of biodiversity. [START_REF] Sampaio | Dinâmica e complexidade da paisagem do Projeto de Assentamento Benfica, Sudeste Paraense[END_REF].

Study on the plant coverage dynamic and on land use, with use of remote sensing, integrating field information allows the definition of parameters that control the spatial and temporal changes which are occurring very quickly in the landscape. Several authors have used this technique to show this process in the Amazon region (including [START_REF] Adams | Classification of multispectral images based on fractions of endmembers: Application to land-cover change in the Brazilian Amazon[END_REF][START_REF] Roberts | ChaAnge identification using multitemporal spectral mixture analysis: Applications in eastern Amazon[END_REF][START_REF] Cochrane | Linear mixture model classification of burned forest in the eastern Amazon[END_REF][START_REF] Souza | An alternative approach for detecting and monitoring selectively logged forests in the Amazon[END_REF][START_REF] Sampaio | Dinâmica na cobertura vegetal e uso da terra do Polígono dos Castanhais no Sudeste Paraense[END_REF]Thalles et al., 2002;[START_REF] Monteiro | Detection of logging in Amazonian transition forests using spectral mixture models[END_REF][START_REF] Souza | Mapping forest degradation in the Eastern Amazon from SPOT 4 through spectral mixture models[END_REF][START_REF] Watrin | Análise evolutiva da cobertura vegetal e do uso da terra em projetos de assentamentos rurais na fronteira agrícola amazônica utilizando geotecnologias[END_REF].

The landscape is considered an applicable scale by researchers from different disciplines to know and assess dynamics of the pioneering fronts in regions with difficult access and fast evolution [START_REF] Venturieri | Utilização de imagens de satélite na caracterização de tipos paisagísticos na frente pioneira de Uruará, Pará[END_REF]. So, the objective is to analyze the expansion of babassu palm (Atalea speciosa Mart. ex Spreng.), in the landscape scale overg eighteen years and its importance for recovery of Legal Reserve Area (LRA) and Permanent Preservation Area (APP) in Benfica SP.

Study area

The area of Benfica-SP is located in the southeastern of Pará State, Itupiranga municipality, between coordinates 0512'20" and 0520'40" South latitude and 4956'40" and 49  48'00" west longitude, at 70 kilometers from the Transamazon Highway and covers an area of 10,026.00 hectares (Figure 1). Benfica II (south of settlement) forming a total of 183 lots occupied since 1994 and have different histories of occupation and uses [START_REF] Dosso | Agriculture ou élevage? Rôle dês couvertures pédologiques dans la différenciation et la transformation de systèmes agraires pionniers au Brésil[END_REF]. Its population is estimated at 1000 inhabitants (Biri Kassoum and Matre d'Hotel, 2002).

184.pdf http://www.isee2012.org/anais/pdf/184.pdf

The original vegetation corresponds to primary forest areas classified according to IBGE (1992) as dense ombrophilous forests and secondary vegetation areas, whose pioneer specie development stages specie vary according to time of abandonment.

Land use consists of acquiring and specializing in cattle raising that causes a progressive substitution of forest by pasture [START_REF] Reynal | Agricultures Familiales et développement em front pionnier amazonien[END_REF] and a period with family agriculture.

Methods

Benfica-SP landscape is a continuous amount of space, segmented in two levels of spatial understanding, perceived in a complimentary and ascendant way: from components (a representative of constituent elements of landscape types) for landscape types, where proportion of areas they occupy is associated to the organization of the community, reflecting its differences and similarities.

The landscape components are compulsory levels in the identification of the plant coverage, image classification, in recognition of the land reality and in the identification of landscape type organization related to description of occupation and use of the land. Its physiognomy informs its function, based on the following criteria: a)

nature, based on its image within the landscape structure, considering its fragmentation; b) proportion (extension), related to representation of the component in terms of the occupied area in the landscape; c) organization in space, according to its shape (regular/ irregular/fragmented) and its distribution (aggregated/disperse, continuous/discontinuous, close/distant), emphasizing the intrinsic organization and interdependence among the components in the system [START_REF] Sampaio | Dinâmica e complexidade da paisagem do Projeto de Assentamento Benfica, Sudeste Paraense[END_REF][START_REF] Laques | Paysage, image et observatoire : lire et diagnostiquer les territoires amazoniens[END_REF].

Landscape type is a combination of components, representative of a portion of the homogeneous space and coherent in relation to the physiognomy of the places (existence of a natural similarity, arrangement and frequency of the constituent elements), both for social and economic use and for ecological operation in the origin of its production.

The study included fieldwork, remote sensing products and techniques, geo- The analysis of dynamic of Benfica-SP landscape types refer to periods : 1987-1992, 1992-1996, 1996-2001 and 2001-2005. TM/Landsat data were processed on the SPRING (Sistema de Processamento de Informações Georreferenciadas) for Windows The most appropriated threshold in the image segmentation values were 8 for similarity and 10 for area (they were selected to define the most suitable limits), whose degree is determined by a parameter t, represented by the Euclidean distance between the vectors associated with each segment. These values were obtained after several approaches until a level of fragmentation considered appropriate for the study.

In classification for areas was used the same limit for analysis of the samples (99.9%), aiming the lowest rejection index, first in the most recent image ( 2005) and then in images of 1987, 1992, 1996 and 2001. In the field, aspects present in the images were correlated with patterns of plant cover and land use. The identified types of correspond to the nature, proportion and organization spatial of its components. The differences in typology established an own configuration, representative of the spatial relationships between different ecosystems or present elements in addition to interactions and natural or anthropic processes (physical, ecological, technical and cultural, etc.).

Representative graphic models of the organization of the landscapes in the field were built, combining components with the productive activities to portray the local and temporal dynamic and to help in the analysis of the landscape changes. Beringuier et al.

(1999) highlighted that these models endeavor to show the spatial organization of the constituent elements of the landscape, so that they are representative of a determined class or typology of the established landscape, by graphic representations.

The delimitation of landscape types in satellite images was based on criteria such as the area of the components, the size of the openings in the forest and their spatial organization.

Results

Landscape components

Ten landscape components were identified (Figure 2 The most studied species in Brazil is Attalea speciosa Mart. ex. Spreng [START_REF] May | Babassu palm in the agroforestry systems in Brazil´s Mid-North region[END_REF], found in isolation, in forests or in open areas, being more frequent in degraded areas, especially in old forest formations cleared in the process of colonization (Ribeiro and Walter, 1988). , 1987-1992, 1992-1996, 1996-2001 and 2001-2005, being Burned 21,54 0,21 0,00 0,00 0,00 0,00 0,00 0,00 0,00 0,00

Landscape dynamics

In dynamics of landscape types (Figure 3), the Forest Landscape (LT0) predominated between 1987 and 1992 with reduction starting in 1996 of almost 80% to 58% in relation to the total area of Benfica-SP.

In the period, the Agricultural Mosaic Landscape (LT1) increased from 6.32% to Palm (LT5) grew by about 60%, 55% and over 100%, respectively.

Discussion

In Southeast of Pará, the expansion of agricultural frontier was followed by pasture implantation starting from periodic deforestation, followed by successive burning. For [START_REF] Santana | Situación y perspectivas de la seguridade alimentaria en la Amazonia: en un marco de cooperacción agropecuária intre-regional[END_REF] first action taken by the settlers is disorderly withdrawal of forest resources for fuel wood, charcoal, poles, stakes and sticks, making it the first source of their subsistence.

In the period 1987-2005, Benfica-SP was marked by many trajectories of evolution that coincided with the productive activities in the process of occupation for pioneer fronts. The arrival of farmer in this period is represented by the absence of capital investment and by a spatial and temporal organization with predominance of Forest Landscape type (TP0).

For [START_REF] Mitja | Renovação das pastagens por agricultores familiares na Amazônia: o caso de Santa Maria, PA[END_REF] these Pioneer fronts consist of a space in formation, where evolutions are very fast and radical. In this context, [START_REF] Piketty | Les déterminants de l'expansion de l'levage bovain en Amazonie Orientale: consquences pour les politiques publiques[END_REF] says that the farmer's productive strategies, based on their culture, traditions and on how they have access to credit and to public policies, were determinant in the organization and dynamic of the landscape.

The great interest of producers by cattle raising in Benfica-SP was established even before their incorporation into land reform. In the period 2001-2005 the expansion of grazing areas with funds from the National Family Agriculture Program (PRONAF), caused environmental and social impacts, with significant reductions in diversity and in areas of forests. Over time, this favored changes in the types of landscapes, from interactions of its components which established ever more or minus, according to conditions and chances of maintaining or not productive systems.

In the context of family agriculture, the expansion of livestock raising system at the cost of using the forest remnants, with time has brought difficulties to maintenance of other productive systems and especially it has brought abandonment and occupation of new areas [START_REF] Silva | Sustentabilidade da agricultura familiar em assentamentos no Sudeste Paraense[END_REF].

Still because of the action of man and his control of pasture, abandonment or possible alternations between the abandonment and use of these managed areas, giving rise to secondary forests of babassu palm (Atalea speciosa Mart. Ex Spreng.) also known as Babassu Palm Forest Landscape (LT4). In this process it is included the financial conditions for labor recruiting.

Babassu palm already present in the forest has spontaneous and abundant regeneration being one of the first palm trees that emerge as the secondary vegetation, mainly along the banks of streams or after withdrawal or burning of original forest, invading damp places, highly degraded, where re-grows vigorously.

The babassu palm tree has an enormous reproductive plasticity that favors its proliferation in pasture areas [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF]. In these areas, fire is the most 184.pdf http://www.isee2012.org/anais/pdf/184.pdf common form of management to eliminate the invaders but it has a contrary effect on the babaçu palm, because after burning the pasture it is precisely the `pindovas´ (young babassu palm individuals) that resist the impact [START_REF] Mitja | Establishment of babassu in pastures in Pará, Brazil[END_REF].

The capitalized producers that have conditions to pay for labor influence land use and pasture remains through management for a longer period preventing faster changes. However, technological conquests present uncertainties from economic point of view, because the high cost of agricultural chemicals elevates cost of recovering the pasture, a fact that often makes the activity of small producers unfeasible [START_REF] Veiga | Expansão e trajetórias da pecuária na Amazônia, Pará, Brasil[END_REF].

In areas of small producers, the babassu palm that normally does not receive crop treatment mixes with the High Fallow (LC4) component and in others the babassu palm spreads dominating old pasture areas. In this process, it forms very often and spontaneously, very dense and dark homogeneous groups, due to the proximity between large palm trees, giving rise to a true Babassu Palm Landscape (LT4).

In turn, in Large Pastures with Babassu Palm Landscape (LT5), where capitalized producers are set down, sometimes the pasture is not abandoned and the babassu palm trees are preserved by proprietor's own will, seeking shade for cattle during the hot hours of the day.

The two examples above is required to comply with environmental laws that determine the use of properties aiming environmental preservation, conservation of natural resources and to ensure a good life quality. Thus, the forests of babassu palm have great importance inside environmental strategies, and provide new economic opportunities as extractivist product.

In the context of environmental management, [START_REF] Homma | Amazônia: como aproveitar os benefícios da destruição?[END_REF] states that the reduction of deforestation and burnings in the Amazon region depends on the development of appropriate economically viable agricultural activities in deforested areas, and the babassu palm is quite suitable for this. According to the Institute of Geography Brazil (IBGE) data, the babassu palm in 2007 was the first in the list of nontimber extractive species contributing significantly to economies of some states in the country.

Thus, Babassu palm has potential for multiple uses and can support the sustainability of family farming systems and the recovery of Legal Reserve Area (LRA)

and Permanent Preservation Areas (PPA) in Benfica SP.

Conclusion

In Benfica-SP, the intensity of use of agricultural systems favored, in time, to strong reactions of Types of Landscape to changes of their components, reflecting different dynamics. This way, recognizing, that both, change is inevitable and the delayed effects that characterize processes in landscape level, the sustainability of the depends on public consistent policies of medium and long term and on the establishment of goals in the management of it.

In the area of Benfica SP, the implantation of agroforestry systems are maintained in the planning, implementation and evaluation of future initiatives certainly would change current scenery of the landscape.

The association of babassu palm with other native species of commercial value rural property would be an alternative to administration of challenges found in pasture management or in the establishment of other cultures, considering the lack of financial resources, and cost reduction and risks for the establishment of exotic species.

The incentive for this activity as a form capable to make changes in the current picture of environment degradation, goes through understanding of lapses found in current productive systems in Benfica-SP and by valorization of a new form of production and handling techniques that meet as much social, economic, political, cultural and ethically, as ecological and environmental principles.

This system could be both made viable in landscapes Babassu Palm However, the establishment of agroforestry systems depends heavily on the commercialization opportunities that vary according with each linked culture. In the case of babassu palm, the use of its full potential, in consortium with other native species as Brazil nuts (Bertholletia excelsa), Copaiba (Copaifera multijuga Hayne), Acai (Euterpe oleracea Mart.) and others could increase food and family incomes generating additional profits.

In addition, the local processing of these species would add value to rural properties and improve social and economic family life. For that, it would be necessary to count on important support policy, mainly for technical assistance, investment 184.pdf http://www.isee2012.org/anais/pdf/184.pdf

Figure 1 .

 1 Figure 1. Location of the area of SP-Benfica, southeastern of Pará State, Itupiranga municipality, Brazil.

  processing and theoretical models. Digital TM-Landsat images were used, orbit/point 224/064, TM3 spectral bands (red), TM4 (close infrared), TM5 (medium infrared) and R5G4B3 colored composition. Images were selected for periods 1987 and 1992 (before occupation of the area), Family Agriculture Program (PRONAF), time of greatest pasture area expansion) and 2005 (before the field surveys started in the second semester of 2005 and first and second semester of 2006).

  version 4.1.1, developed by the Brazilian Space Research Institute (INPE). The classification method supervised by regions associated the image resulting from the process of image segmentation for fragmentation into homogeneous units.

  The components Pasture with Woody Species (LC7), Pasture with Babassu Palm (LC8) and Clean Pasture (LC10) consist mainly of species such as Brachiaria brizantha (Hochst. ex A. Rich.) Stapf. and Panicum sp. The Clean Pasture corresponds to recently implanted or well-managed areas, with a low degree of invaders; the Pasture with Woody Species presents invaders of various species, but still sustains pasture activity and Pasture with Babassu Palm is dominated by young Babassu palm individuals. Those components were observed for 18 years in four different periods

Figure 2 .

 2 Figure 2. Landscape components of the SP-Benfica, southeastern of Pará State, Itupiranga municipality, Brazil.

  There have been identified six landscape types as: Forest Landscape (LT0), Agricultural Mosaic Landscape (LT1), Agricultural Mosaic with Pasture (LT2) Large Pastures Landscape (LT3), Babassu Palm Forest Landscape (LT4) and Large Pastures with Babassu Palm (LT5). In the landscape (LT0), the Forest Remnant component (LC1) occupies the landscape in an almost continuous way and in the Agricultural Mosaic (LT1) very small portions of land are organized in a heterogeneous manner, generally rice, corn and cassava fields or rice with pasture consortium. These practices occur shortly after felling and burning the forest, with the surroundings areas occupied by High Fallow (LC4) and Low Fallow (LC6) components. The Agricultural Mosaic with Pasture Landscape (LT2) presents a large area of forest close to anthropic areas, diversity of components and the agricultural production systems are giving place to areas occupied by fallow because of abandonment of fields and pasture. The Large Pasture Landscape (LT3) corresponds to a matrix of pasture consisted of Clean Pasture (LC10) and/or Pasture with Woody Species (LC7), where direct 184.pdf http://www.isee2012.org/anais/pdf/184.pdf substitution of forest by pasture is executed. The isolated and distant remnants are substituted by new pasture, during a new agricultural cycle. When the cover is not totally formed by pasture, there is dominance of Low Fallow (LC6) component. Signaling management of the area, rarely is observed the presence of adult and young babassu palm trees or taboa (Typha sp.), a species that prefers wet places such as swamp and whose presence becomes problematic for farmers for being very competitive. The Babassu Palm Forest Landscape (LT4) is characterized by a matrix where the Babassu Palm Forest (LC9) and High Fallow (LC4) components predominate, indicating abandonment of old pasture, and the Large Pasture with Babassu Palm Landscape (LT5), reflects the expansion of the component Pasture with Babassu Palm (LC8).

  14.79% giving room for Agricultural Mosaic with Pasture Landscape (LT2) that increased from 11.84% to 23.99%. The course of this landscape (LT2) to the Large Pasture Landscape (LT3) in the 1987-1992 period presented a reduction in the area from 2.15% to 1.54%, and may be associated to abandonment of management in these areas and to increase in the area of Agricultural Mosaic with Pasture Landscape (LT2).From 1992 to 1996 the area of Forest Landscape (LT0) reduced 19% contributing to the expansion of the Large Pasture Landscape (LT3) and Large Pasture with Babassu Palm Landscape (LT5) and to the appearance of Babassu Palm Landscape (LT4).In period 1996-2001, the largest changes were registered in landscape types. In addition to disappearance of the Forest Landscape (LT0), the Agricultural Mosaic Landscape (LT1), the Agricultural Mosaic with Pasture Landscape (LT2) and Large Pasture Landscape (LT3) increased by over 80%, 70% and 174%, respectively. In this period the Babassu Palm Forest Landscape (LT4) increased by almost 20% and the Babassu Palm component (LC9) dominated more than 65% of its area.

  Figure 3. Dynamics of landscape types in SP-Benfica, southeastern of Pará State, Itupiranga municipality, Brazil.

  , Piketty et al. (2005), Poccard-Chapuis et al. (2005) e Sartre et al. (2005) the importance of livestock is unquestionable in the regional economy. In this context,[START_REF] Silva | Manejo Florestal Sustentável da Caatinga: Adequação Ambiental e Produção de Energia Agroecológica em Projetos de Assentamento e Propriedades Coletivas do Plano Nacional de Reforma Agrária[END_REF] highlights that, traditionally, the 184.pdf http://www.isee2012.org/anais/pdf/184.pdf

  Forest(LT4) and Large Pastures with Babassu Palm (LT5) formed by abandonment of pasture where predominates the natural regeneration of babassu palm, as in landscape Large Extensions Pasture (LT3) where the babassu palm is suppressed during pasture management.

  Typha sp., Ludwigia latifolia (Benth.) H. Hara and Panicum pilosum Sw.

	): Forest Remnant (LC1)
	referring to the remains of the original forest, that still present a great species diversity;
	Riparian Forest (LC2), forest areas on the banks of bodies of water, sometimes well
	conserved and other times disturbed; Swamp (LC3) secondary vegetation derived from
	degradation of Riparian Forests and presenting a great amount of secondary species
	such as Secondary forests present different successional stages, that depend on time of
	abandonment and on the history of use of the area, as: High Fallow (LC4) refers to the 8

to 15 year-old secondary forests, approximately; Fallow with Jurubeba (LC5), young fallow, average age 5 years, dominated by Solanum rugosum Dunal, commonly called Jurubeba; Low Fallow (LC6), young fallow, average age five years, without apparent dominance by any species and Babassu Palm (LC9), secondary forests dominated by babassu palm with age around 20 years old. Babassu is the common name for several species of the genus Attalea (A.speciosa A. brasiliensis), pionner species of Amazon forest and of Cerrado biome.

TABLE 1 .

 1 Quantification of areas the Landscape Components of the SP-Benfica, southeastern of Pará State, Itupiranga municipality, Brazil for years 1987,1992, 1996, 2001 and 2005.

	Components	1987	1992		1996		2001		2005
	Landscape	ha	%	ha	%	ha	%	ha	%	ha	%
	(LC1)	8605,4	85,83 7810,87 77,91 6709,19 66,92 3772,10 37,62 2486,62 24,80
	(LC2)	3 211,59 2,11 200,86 2,00 326,97 3,26 134,61 1,34	83,72	0,84
	LC3)	36,17 0,36	14,56	0,15	49,60	0,49	21,80	0,22	98,08	0,98
	LC4)	321,89 3,21 520,90 5,20 788,08 7,86 1586,04 15,82 1009,81 10,07
	(LC5)	1,76	0,02 149,82 1,49	60,02	0,60 441,55 4,40 519,77 5,18
	(LC6)	193,96 1,93 323,03 3,22 543,36 5,42 1219,25 12,16 1465,38 14,62
	(LC7)	217,90 2,17 137,18 1,37 212,24 2,12 704,05 7,02 1536,23 15,32
	(LC8)	55,53 0,55 383,23 3,82 484,55 4,83 786,73 7,85 1165,62 11,63
	(LC9)	36,81 0,37 157,56 1,57 427,23 4,26 536,31 5,35 831,53 8,29
	(LC10)	265,00 2,64 326,02 3,25 306,98 3,06 606,99 6,05 775,49 7,73
	Other: Water 0,00	0,00	0,00	0,00	26,96	0,87	6,91	0,07	3,68	0,04
	Exposed Soil 58,43 0,58	1,97	0,02	90,82	0,31 209,65 2,09	50,06	0,50

02/04/13 00:51

sur 17 02/04/13 00:51

184.pdf http://www.isee2012.org/anais/pdf/184.pdf 4 sur 17 02/04/13 00:51

sur 17 02/04/13 00:51

184.pdf http://www.isee2012.org/anais/pdf/184.pdf

Acknowledgements

This project has benefited from the support of IFB (Institut Français de la Biodiversité), of IRD (Institut de Recherche pour le Développement), of CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico), from the UFRA (Universidade Federal Rural da Amazônia). The authors would like to deeply thank the farmers of Benfica-SP community for their help, availability and conviviality during the accomplishment of this work.

Code from now on recovery.

Although the New Forest Code establish restrictions on use of Legal Reserve (shallow cutting, land use alteration and exploitation for commercial purposes, with some exceptions), it makes possible the exploration through sustainable management plans.

credits, processing and commercialization, technology research and industrialization of production beyond expansion of infrastructure and organization of farmers.

In ecological and environmental context, sustainable use of babassu in agroforestry systems could contribute to recovery actions of the Legal Reserve Areas (LRA) and Permanent Preservation Areas (PPA), which are a priority for the mitigation of environmental liabilities and thus, would meet the requirements of the New Forest