Jawhar Chebbi

Vincent Dubanchet

José Alvaro Perez Gonzalez

Daniel Alazard
email: alazard@isae.fr

Linear dynamics of flexible multibody systems A system-based approach

Keywords: Flexible structure, Multibody system, Linear system B D. Alazard

We present a new methodology to derive a linear model of flexible multibody system dynamics. This approach is based on the two-port model of each body allowing the model of the whole system to be built just connecting the inputs/outputs of each body model. Boundary conditions of each body can be taken into account through inversion of some input-output channels of its two-port model. This approach is extended here to treat the case of closed-loop kinematic mechanisms. Lagrange multipliers are commonly used in an augmented differential-algebraic equation to solve loop-closure constraints. Instead, they are considered here as a model output that is connected to the adjoining body model through a feedback. After a summary of main results in the general case, the case of planar mechanisms with multiple uniform beams is considered, and the two-port model of the Euler-Bernoulli beam is derived. The choice of the assumed modes is then discussed regarding the accuracy of the first natural frequencies for various boundary conditions. The overall modeling approach is then applied to the well-known four-bar mechanism.

Introduction

Dynamics of flexible multibody systems has been vigorously studied over the last decades in connection with various application fields such as light weight robotic arm, rotating machinery, helicopter rotor with flexible blades, large flexible space structures, and so on. We can find hundreds of papers about the writing, solving, and simulation of the differential and algebraic equations governing the dynamics of flexible multibody systems [START_REF] Schiehlen | Multibody system dynamics: roots and perspectives[END_REF][START_REF] Shabana | Flexible multibody dynamics: review of past and recent developments[END_REF]. We can distinguish methods used to build a validation model including nonlinear terms (centrifugal, Coriolis) [START_REF] De Luca | Closed-form dynamic model of planar multilink lightweight robots[END_REF][START_REF] Mohan | A recursive, numerically stable, and efficient simulation algorithm for serial robots with flexible links[END_REF] and methods used to derive a preliminary design model, which can be also used for control design [START_REF] Krauss | Transfer matrix modeling of systems with noncollocated feedback[END_REF][START_REF] Murali | Mechanical-attitude controller co-design of large flexible space structures[END_REF]. The latter often focus on the linear behavior of flexible multibody systems under a small displacement assumption. Link flexibility is approximated by discretized models using finite element methods (FEM) or assumed mode methods (AMM) [START_REF] Theodore | Comparison of the assumed modes and finite element models for flexible multilink manipulators[END_REF]. The model of the whole multibody system is then derived using the Euler-Lagrange formalism. For example, in space applications, considering a flexible robotic arm fitted on a chaser spacecraft, the inertial load seen from the tip of a flexible link can vary considerably according to the geometrical configuration of the robotic arm. The AMM is then often used [START_REF] De Luca | Closed-form dynamic model of planar multilink lightweight robots[END_REF][START_REF] Theodore | Comparison of the assumed modes and finite element models for flexible multilink manipulators[END_REF] with the first clamped-mass eigenfunctions. The main drawback is that the mass must be a time-variant parameter and need to be updated along with the geometrical configuration. Otherwise, an approximation by an average value is still possible. However, if we consider the capture of a massive space debris, then the boundary conditions at the tips of each flexible link may be closer to the clamped-clamped condition than to the clamped-mass condition. During the preliminary design study of such a system, we can argue that the system dynamics is more sensitive to such a change in the boundary conditions than in the effects of nonlinear terms. Thus, there is a real need to have a model for each link or substructure that is valid for arbitrary boundary conditions before doing the assembly of the whole structure.

Among all substructuring approaches, the transfer matrix method (TMM) has motivated lots of literature during last decades. Leckie [START_REF] Leckie | Transfer-matrix fundamentals[END_REF] stated the fundamentals of this method: the transfer matrix is a relation between the left point and right point state vectors of a flexible body. The state vector at a point is defined as the augmented vector composed of the generalized displacements and the corresponding generalized forces. This method was particularly well suited for modeling structures composed of serially connected bodies or open-chainlike structures. Later, this approach was merged with the finite-element (FE) method and presented as the FE-TM method [START_REF] Dokainish | A new approach for plate vibrations: combination of transfer matrix and finite-element technique[END_REF]. The FE-TM method was first developed in [START_REF] Mucino | An exact condensation procedure for chain-like structures using a finite element-transfer matrix approach[END_REF] to provide reduced models using a condensation procedure, before being extended in [START_REF] Rong | Modified finite element transfer matrix method for eigenvalue problem of flexible structures[END_REF] for complex shape bodies to reduce computation efforts in solving eigenvalue problems. Some extensions were also developed for control design purposes [START_REF] Tan | A modified finite element-transfer matrix for control design of space structures[END_REF], including noncollocated feedback [START_REF] Krauss | Transfer matrix modeling of systems with noncollocated feedback[END_REF], and for hybrid (rigid and flexible) multibody systems [START_REF] Rui | Transfer matrix method for linear multibody system[END_REF]. In [START_REF] Kitis | Natural frequencies and mode shapes of flexible mechanisms by a transfer matrix method[END_REF], the TMM was applied to a closed-loop kinematic chain of flexible bodies and more particularly to the four-bar mechanism. Nevertheless, the TMM has some limitations due to the inversion of submatrices, which are not always square or invertible, depending on the boundary conditions [START_REF] Tan | A modified finite element-transfer matrix for control design of space structures[END_REF]. Furthermore, another limitation comes from the fact that this method cannot be directly applied to tree-like structures, like a spacecraft composed of a main rigid body fitted with several flexible appendages (solar panels, antenna, etc.). In such a case, the objective of the modeling is to find a relation between the generalized forces and displacements at the root of the tree-like structure, that is, the main body. Such a model is required for the design of the attitude and orbit control system (AOCS). Some alternatives and efficient approaches were developed to treat that case. They are based on the effective mass/inertia model of each flexible appendage (i.e., the relation between forces and accelerations at the anchorage point of the appendage on the main body) [START_REF] Girard | Dynamics of Structure in Industry[END_REF] or also on the impedance matrix [START_REF] Pascal | Dynamics analysis of a system of hinge-connected flexible bodies[END_REF][START_REF] Sylla | Dynamics of a rotating flexible and symmetric spacecraft using impedance matrix in terms of the flexible appendages cantilever modes[END_REF]. These methods were first introduced in [START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF] to analyze dynamic coupling between substructures through the component mode synthesis (CMS) and are commonly used in space engineering [START_REF] Guy | Dynamic modeling and analysis of spacecraft with variable tilt of flexible appendages[END_REF]. The limitation of the effective mass/inertia approach is the loss of appendage information beyond the anchorage point, such as the deflection at the free tip. Therefore, such a method cannot be applied to model an arbitrary chain-like structure of flexible bodies.

The main contribution of this paper is to propose a system-based approach to model the linear dynamic behavior of a flexible link in a multibody system. This approach is independent of the boundary conditions applied at the connecting points of this link. The inputs and outputs of this model, named two-input-two-output port (TITOP) model, are the forces Fig. 1 Multibody systems and accelerations at the connecting points of the link with its neighbors. This model embeds, in the same minimal state-space representation, the direct and the inverse dynamic models of the link. Its various input-output channels (or ports) are thus invertible. By applying very simple channel inversion operations this model can be directly plugged into the whole multibody system model, seen as a block-diagram model. The boundary conditions at a connecting point of the link are seen as an external feedback loop between forces and accelerations. The same model of the link can be used for any kind of open-loop or closedloop kinematic mechanisms involving this link. This method fills the gap between the TM method and the effective mass-inertia method. Note that from the terminology point of view, what is called the state vector of left and right points in the TM approach are in fact the inputs and outputs of the transfer matrix. These inputs and outputs are completely revisited in the TITOP approach. Finally, accelerations at the connecting points are considered instead of generalized displacements. The state vector is now the vector of variables associated to the minimal state-space representation of the TITOP model (which is twice the number of flexible modes of the considered body). First introduced in [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF], the TITOP model approach is also illustrated on space engineering applications in [START_REF] Murali | Mechanical-attitude controller co-design of large flexible space structures[END_REF][START_REF] Perez | Linear dynamic modeling of spacecraft with open-chain assembly of flexible bodies for ACS/structure co-design[END_REF][START_REF] Perez | A flexible appendage model for use in integrated control/structure spacecraft design[END_REF][START_REF] Perez | Flexible multibody system linear modeling for control using component modes synthesis and double-port approach[END_REF].

In Sect. 2, the main results of [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF] are summarized and reformulated in the context of general multibody system with a focus on closed-loop mechanisms. In Sect. 3, the TITOP model is derived for an Euler-Bernoulli beam. It will be shown that the first two natural frequencies for any kind of tip conditions are accurately modeled from the simple TITOP approach. In Sect. 4, the methodology is then applied to the well-known four-bar mechanism: the first natural frequencies and modal shapes are analyzed versus the crank angle and compared with the results already published [START_REF] Kitis | Natural frequencies and mode shapes of flexible mechanisms by a transfer matrix method[END_REF], [START_REF] Turcic | Dynamic analysis of elastic mechanism systems. Part I: applications[END_REF].

Two-port model approach

An N degrees-of-freedom (d.o.f.) multibody system can be divided into bodies (or links) connected to each other by (prismatic, revolute, cardan, spherical) joints. We can distinguish chain-like multibody systems (Fig. 1a), tree-like systems (Fig. 1b), and closed-loop kinematic systems (Fig. 1c). In all the cases, the objective of this paper is to derive the linear model of the system dynamics around a given equilibrium configuration of N d.o.f., Fig. 2 The ith flexible link in the structure θ = [θ 1 , θ 2 , . . . , θ N ] T and θ = 0, assuming small variations of the d.o.f. around this configuration and small deformations of flexible links. Therefore, the Coriolis and centrifugal forces are neglected throughout this paper.

Frame definition Under this small displacement assumption, it is possible to define an inertial frame R i = (P 0 i , x i , y i , z i ) attached to the equilibrium configuration of each link i. In the chain-like and tree-like cases, P 0 i is the connection point of link i with the parent link at the equilibrium. P i/j is the direction cosine matrix (DCM) between frames R i and R j (i.e., the matrix of components of unitary vectors x j , y j , z j in R i ). P i/j depends only on the given configuration θ for all i and j .

Two-port model of a link

The link L i connected to the parent substructure L i-1 at the point P i and to the child substructure L i+1 at point C i is depicted in Fig. 2. The double-port or TITOP (two-input-twooutput port) model G L i P i ,C i (s), proposed in [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF], is a linear dynamic model between 12 inputs: -the six components in R i of the wrench

W C i = [F T C i T T C i ]
T applied by the substructure L i+1 to the link L i at point C i : F C i stands for the three-component force vector applied at point C i , and T C i stands for the three-component torque vector applied at point C i , -the six components in R i of the acceleration twist ẍP i = [a T P i ωT P i ] T (time-derivative of the twist) of point P i : a P i stands for the three-component linear acceleration vector at point P i , and ωP i stands for the three-component angular acceleration vector at point P i , and 12 outputs:

-the six components in R i of the acceleration twist

ẍC i = [a T C i ωT C i ] T of point C i , -the six components in R i of the wrench W P i = [F T P i T T P i ]
T applied by the link L i to the substructure L i-1 at point P i , and can be represented by the block-diagram depicted in Fig. 3.

The way to obtain such a TITOP model G L i P i ,C i (s) will be detailed in Sect. 3 for Euler-Bernoulli beams and is also detailed in [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF] and [START_REF] Perez | Flexible multibody system linear modeling for control using component modes synthesis and double-port approach[END_REF] in the general case. The general procedure is based on the modal analysis of the link L i with the clamped at P i -free at C i boundary

L i P i ,C i (s) of link L i
conditions. Indeed, the Lagrange formulation of dynamics using the assumed mode method (AMM) or finite element method (FEM) leads to a generalized second-order differential equation [START_REF] Shabana | Flexible multibody dynamics: review of past and recent developments[END_REF][START_REF] Theodore | Comparison of the assumed modes and finite element models for flexible multilink manipulators[END_REF]:

D L i P i L T P i L P i 1 n i ẍP i ηi + 0 0 0 diag(2ξ j ω j ) ηi + 0 0 0 diag(ω 2 j ) η i = -I 6 τ T C i P i 0 Φ T C i W P i W C i , (1)
where:

n i is the number of flexible clamped-free modes of link L i characterized by the modal coordinates vector η i , the frequencies ω j , j = 1, . . . , n i , and the damping ratio ξ j , j = 1, . . . , n i , -1 n i is the identity matrix of size n i , -L P i is the n i × 6 modal participation factor matrix of the link at point P i and projected in the frame R i , -Φ C i is the 6 × n i projection matrix of the n i clamped-free modal shapes on the six d.o.f. at point C i and projected in the frame R i , τ C i P i is the "rigid" kinematic model between point C i and P i projected in the frame R i :

τ C i P i = 1 3 ( * --→ C i P i ) 0 3×3 1 3
, where ( * --→ C i P i ) is the skew-symmetric matrix associated with the vector from C i to P i , -D L i P i is the 6 × 6 rigid mass model of the link at point P i and projected in the frame R i . The right-hand term of Eq. (1) describes the contribution of the wrenches W P i and W C i to the generalized force vector. The acceleration twist ẍC i at point C i can be easily expressed from the generalized acceleration vector, the modal shapes, and the kinematic model: and:

ẍC i = [τ C i P i Φ C i ] ẍP i ηi (2)
From Eqs. ( 1) and ( 2), the general state-space realization (see also the definition of a state-space realization of a given transfer in Appendix 1) of the TITOP model G

L i P i ,C i (s) reads: ⎡ ⎢ ⎢ ⎣ ηi ηi ẍC i W P i ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 0 n i ×n i 1 n i 0 n i ×6 0 n i ×6 -diag(ω 2 j ) -diag(2ξ j ω j ) Φ T C i -L P i -Φ C i diag(ω 2 j ) -Φ C i diag(2ξ j ω j ) Φ C i Φ T C i (τ C i P i -Φ C i L P i ) L T P i diag(ω 2 j ) L T P i diag(2ξ i ω i ) (τ C i P i -Φ C i L P i ) T -D L i P i + L T P i L P i ⎤ ⎥ ⎥ ⎥ ⎦ × ⎡ ⎢ ⎢ ⎣ η i ηi W C i ẍP i ⎤ ⎥ ⎥ ⎦ . ( 3
)
This model embeds, in the same minimal state-space realization, the direct dynamic model (transfer from acceleration to force) of link i at the point P i and the inverse dynamic model (transfer from force to acceleration) of link i at the point C i . All data required in Eq. ( 3) are commonly provided by FEM software. In [START_REF] Perez | Flexible multibody system linear modeling for control using component modes synthesis and double-port approach[END_REF], the link between the TITOP model and substructuring methods like the Graig-Bampton decomposition and the component modes synthesis (CMS) is detailed in the general case. In [START_REF] Murali | Mechanical-attitude controller co-design of large flexible space structures[END_REF], a custom NASTRAN/SIMULINK interface was developed to declare and manipulate the TITOP model G L i P i ,C i (s) as a linear dynamic system object under the MATLAB/SIMULINK environment.

Let us note that the modal analysis of the left-hand term of Eq. ( 1) corresponds to the free (at P i )-free (at C i ) case and provides six rigid modes in addition to the free-free flexible modes. These rigid modes are not taken into account in the model G L i P i ,C i (s) defined by Eq. ( 3). This is one the main advantages of the TITOP model: the boundary conditions and the rigid modes can be taken into account "outside" the model. This is detailed in the following section.

Interconnection and channel inversion of TITOP models

The open-loop dynamics of the TITOP model G L i P i ,C i (s), that is, when the inputs are null, corresponds to the link dynamics when it is clamped at P i (ẍ P i = 0) and free at C i (W C i = 0). It must be noticed that all the 12 input-output channels of the model G L i P i ,C i (s) are invertible assuming that a residual mass (inertia) of link i is attached to the points P i and C i along (around) the three axes. Such an assumption is valid by using the finite-element method or the assumed mode method to model the flexibility of the link.

[G

L i P i ,C i ] -1 I (s)
denotes the model where the channels numbered in the index vector I are inverted according to the procedure described in Appendix 1. This channel inversion operation is very useful to analyze the dynamics of the link L i for various boundary conditions:

-clamped at P i and clamped at C i , using [G [1:6] (s) whose inputs are ẍC i and ẍP i , -free at P i and free at C i , using [G L i P i ,C i ] -1 [7:12] (s) whose inputs are W C i and W P i , -free at P i and clamped at C i , using [G L i P i ,C i ] -1 (s) whose inputs are ẍC i and W P i , -or to take into account a revolute (resp. prismatic) joint at the point P i around (resp. along) the z i -axis using [G

L i P i ,C i ] -1
L i P i ,C i ] -1 12 (s) (resp. [G L i P i ,C i ] -1 9 (s))
, that is, no torque (resp. force) can be applied by the link L i to the parent substructure L i-1 around (resp. along) the z i -axis. A spherical joint can be taken into account using [G [10:12] (s). More generally, once the TITOP model of each link is defined and the DCM P i/j for a given nominal configuration θ are computed, the linear model of the whole multibody system can be built by simply connecting the inputs/outputs of the TITOP models. Hereafter, some examples of interconnection are presented.

L i P i ,C i ] -1

Chain-like multibody system

Let us consider a flexible spacecraft S whose reference body frame is (O, x 0 , y 0 , z 0 ), fitted with a flexible boom B at point P 1 with a revolute joint along the (P 1 , z 1 )-axis. The boom holds a flexible antenna A at point P 2 . A sketch of such a multibody system is depicted in Fig. 4 (left). The block-diagram associated to this system is presented in Fig. 4 (right). It involves the three TITOP models of the three bodies with the channel inversions required by the different boundary conditions. P a i/j is the augmented 6 × 6 DCM: In this block-diagram, seven double integrators are added on the outputs to represent the seven rigid modes associated to the seven d.o.f. of this system: the six d.o.f. of the spacecraft and one d.o.f. of the revolute joint. The seven inputs of this model are the torque u 1 applied inside the revolute joint by the driven mechanism and the six components of the wrench W 0 applied by the attitude and orbit control system at the reference point O. All dynamic couplings between the flexible modes of the antenna, the boom, and the spacecraft are directly taken into account through the feedback connections between the three TITOP models. Furthermore, the inertia j 1 , the stiffness k 1 , and the damping factor f 1 of the driven mechanism located in the revolute joint can be added by a feedback loop u 1 = -j 1 δθ 1f 1 δθ 1k 1 δθ 1 , directly plugged into the block-diagram. This can be used to handle any arbitrary boundary condition in the revolute joint, from the free condition (

P a i/j = P i/j 0 3×3 0 3×3 P i/j .
j 1 = k 1 = f 1 = 0) to the clamped condition (k 1 → ∞).
Tree-like multibody system In this case, one link B may have two (or more) child bodies connected to B at points C 1 and C 2 . This case is quite similar to the previous one and requires a model of the link augmented with additional ports. This case is not detailed in this paper. The reader is advised to refer to [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF] where the three-port model of a link (denoted

G B P ,C 1 ,C 2 (s))
and an example of a tree-like system are detailed (see Fig. 18 in [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF]).

Closed-loop multibody system Such mechanisms are commonly modeled through the Lagrange approach using Lagrange multipliers [START_REF] Shabana | Flexible multibody dynamics: review of past and recent developments[END_REF][START_REF] Simeon | On Lagrange multipliers in flexible multibody dynamics[END_REF]. They are used to be denoted λ and correspond to the generalized constraint forces. Considering the model of the link L i given in Eqs. ( 1) and ( 2) with a kinematic-chain closed at point C i , the Lagrange multipliers are nothing else than the wrench at this point applied by the link L i : λ = -W C i . The classical Lagrange approach leads to an augmented system of differential-algebraic equations: Considering the TITOP formulation of the model in Eq. ( 3), λ (up to the sign) is an input of G

⎡ ⎣ D L i P i L T P i τ T C i P i L P i 1 n i Φ T C i τ C i P i Φ C i 0 ⎤ ⎦ ⎡ ⎣ ẍP i ηi λ ⎤ ⎦ + ⎡ ⎣ 0 0 0 0 diag(2ξ j ω j ) 0 0 0 0 ⎤ ⎦ ⎡ ⎣ ηi ⎤ ⎦ + ⎡ ⎣ 0 0 0 0 diag(ω 2 j ) 0 0 0 0 ⎤ ⎦ ⎡ ⎣ η i ⎤ ⎦ = ⎡ ⎣ -W P i 0 ẍC i ⎤ ⎦ . ( 4
)
L i P i ,C i (s) or an output of the model [G L i P i ,C i ] -1 [1:6] (s)
where the upper port is inverted. Thus, the upper port of [G [1:6] (s) can be connected in a feedback loop to the upper port of the adjoining link in order to close the kinematic chain. Such an interconnection is presented in Fig. 5 (right) in the case of a very simple closed-loop structure composed of two links L 1 and L 2 clamped to the ground and one to each other at point C 1 ≡ C 2 . In that case, the model outputs are the two wrenches W P 1 and W P 2 applied by the system to the ground at points P 1 and P 2 . Similarly to the case of the chain-like multibody system example, a revolute joint around the (C 2 , z 1 )-axis can be taken into account using the model [G [1:5] (s) to cancel the torque around this axis in the joint between L 1 and L 2 , that is, λ(6) = W C i (6) = 0 (see also the illustration on the four-bar mechanism in Sect. 4).

L i P i ,C i ] -1
L 2 P 2 ,C 2 ] -1

Two-port model of a beam in the planar case

Let us consider a uniform beam B between points P and C, characterized by the following parameters (see Fig. 6):

-mass density ρ (kg/m 3 ), -section S (m 2 ), -length l (m), -Young modulus E (N/m 2 ), -second moment of area w.r.t. the z-axis: I z (m 4 ).

The objective is to compute the TITOP model G B P ,C (s). In the (3 d.o.f.) planar case, this model will be denoted G(s) and is a transfer between six inputs:

-three components F x C (N), F y C (N), and T z C (Nm) in the frame R = (P 0 , x, y, z) of the planar wrench applied to the beam at point C, -three components üC (m/s 2 ), vC (m/s 2 ), and αC (rad/s 2 ) in the frame R of the planar acceleration twist of the beam at point C, -three components F x P (N), F y P (N), and T z P (Nm) in the frame R of the planar wrench applied by the beam at point P . This model G(s) is derived in following sections considering a decoupling between the pure flexion in the plane (P 0 , x, y) and the pure traction-compression along the (P 0 , x)axis.

Pure flexion in the plane (P 0 , x, y)

A well-known approach to derive the dynamic model of a flexible beam is based on the Lagrange technique combined with the assumed modes method (AMM) [START_REF] De Luca | Closed-form dynamic model of planar multilink lightweight robots[END_REF]. This approach involves a moving body frame R = (P , x, ỹ, z) attached to the beam at point P (see Fig. 7). Under the small deflection assumption, the displacement v(x, t) along the y-axis at any point of abscissa x and at any time t reads as follows:

v(x, t) = v P (t) + xα P (t) + ṽ(x, t) ∀x ∈ [0, l].
(

) 5
The deflection ṽ(x, t) in the moving frame R is split up into n shape functions

φ i (x) (φ(x) = [φ 1 (x) • • • φ n (x)] T being the shape function vector), associated with n time-domain functions q i (t), n = 1, . . . , n (q(t) = [q 1 (t) • • • q n (t)] T ): ṽ(x, t) = n i=0 φ i (x)q i (t) = φ(x) T q(t). (6)
Then, the kinetic energy reads

T = 1 2 ρS l 0 vP + x αP + Σ n i=0 φ i (x) qi 2 dx = 1 2 ⎡ ⎣ vP αP q ⎤ ⎦ T M rr M T rf M rf M ff ⎡ ⎣ vP αP q ⎤ ⎦ ,
where:

- The elastic potential energy reads

M rr = ρSl ρSl 2
V = 1 2 EI z l 0 ∂ 2 ṽ(x, t) ∂x 2 2 dx = 1 2 ⎡ ⎣ v P α P q ⎤ ⎦ T 0 2×2 0 2×n 0 n×2 K n×n ⎡ ⎣ v P α P q ⎤ ⎦ ,
where

K = EI z l 0 [ d 2 φ(x) dx 2 ][ d 2 φ(x) dx 2 ] T dx = EI z l 0 φ (x)φ T (x) dx is the n × n stiffness ma- trix.
Applying the action-reaction principle at point P , the work of external loads is

W ext = -v P F y P -α P T z P + v(l)F y C + v (l)T z C (7) = -v P -α P v P + lα P + q T φ(l) α P + q T φ (l) ⎡ ⎢ ⎢ ⎣ F y P T z P F y C T z C ⎤ ⎥ ⎥ ⎦ (8)
= ⎡ ⎣ v P α P q ⎤ ⎦ T ⎡ ⎣ -1 0 1 0 0 -1 l 1 0 n×1 0 n×1 φ(l) φ (l) ⎤ ⎦ ⎡ ⎢ ⎢ ⎣ F y P T z P F y C T z C ⎤ ⎥ ⎥ ⎦ , ( 9
)
and the Lagrange derivation leads to the following second-order equation:

M rr M T rf M rf M ff ⎡ ⎣ vP αP q ⎤ ⎦ + 0 2×2 0 2×n 0 n×2 K n×n ⎡ ⎣ v P α P q ⎤ ⎦ = -I 2 τ T 0 n×2 Φ T (l) ⎡ ⎢ ⎢ ⎣ F y P T z P F y C T z C ⎤ ⎥ ⎥ ⎦ (
From Eqs. [START_REF] Krauss | Transfer matrix modeling of systems with noncollocated feedback[END_REF] and [START_REF] Leckie | Transfer-matrix fundamentals[END_REF] we can easily derive the state-space representation, associated with the state vector x f = [q T qT ] T , of the TITOP model G f (s) relative to the pure flexion dynamics in the plane (P 0 , x, y) (see Fig. 8):

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ q q vC αC F y P T z P ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 0 n 1 n 0 n×2 0 n×2 -M -1 ff K 0 n M -1 ff Φ T (l) -M -1 ff M rf -Φ(l)M -1 ff K 0 n Φ(l)M -1 ff Φ T (l) τ -Φ(l)M -1 ff M rf M T rf M -1 ff K 0 n (τ -Φ(l)M -1 ff M rf ) T M T rf M -1 ff M rf -M rr ⎤ ⎥ ⎥ ⎥ ⎦ × ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ q q F y C T z C vP αP ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 12
)

Choice of the modal decomposition

In this section, we discuss the choice of the decomposition defined in Eq. ( 6) regarding the capability of the TITOP model G f (s) (and its inverses) to take into account various boundary conditions. Among all possible boundary conditions, a very simple test consists in computing and comparing the natural frequencies ω cf,i and ω f c,i of models G f (s) and G -1 f (s), respectively. Indeed, G f (s) is the model of the beam under the clamped-free condition and G -1 f (s) is the model of the beam under the free-clamped condition. Since the beam is uniform, the natural frequencies of the two models must be equal: ω cf,i = ω f c,i ∀i.

Two decompositions are considered to build the TITOP model G f (s):

-the well-known clamped-free assumed mode decomposition [START_REF] Mohan | A recursive, numerically stable, and efficient simulation algorithm for serial robots with flexible links[END_REF]. The model obtained from this decomposition is denoted G f,n (s), where n is the number of assumed modes. -the decomposition based on fifth-order polynomial shape functions [START_REF] Murali | Mechanical-attitude controller co-design of large flexible space structures[END_REF]. The model obtained from this decomposition is denoted G f,pol (s).

Model G f,n (s)

In this case, the decomposition is based on the n first modes of the Euler-Bernoulli beam in the clamped-free condition, and the shape functions in (6) are [START_REF] Bishop | The Mechanics of Vibration[END_REF][START_REF]Euler-Bernoulli beam theory[END_REF]

φ i (x) = cos β i x -cosh β i x - cos β i l + cosh β i l sin β i l + sinh β i l (sin β i x -sinh β i x), ( 13
)
where

β i , i = 1, . . . ,

Model G f,pol (s)

In this approach, given in more details in [START_REF] Murali | Mechanical-attitude controller co-design of large flexible space structures[END_REF], the deflection ṽ(x, t) in the moving frame R is expressed as ṽ(x, t) = a 2 (t)x 2 + a 3 (t)x 3 + a 4 (t)x 4 + a 5 (t)x 5 . This fifth-order polynomial ensures that ṽ(0, t) = ṽ (0, t) = 0 ∀t . The coefficients a i (t) are then expressed as a function of four time-domain functions q i (t) defined by q(t) = ṽ (0, t) ṽ(l, t) ṽ (l, t) ṽ (l, t) T , [START_REF] Murali | Mechanical-attitude controller co-design of large flexible space structures[END_REF] and the shape functions vector reads

φ(x) = ⎡ ⎢ ⎢ ⎣ 1 2 x 2 -3 2l x 3 + 3 2l 2 x 4 -1 2l 3 x 5 10
l 3 x 3 -15 l 4 x 4 + 6 l 5 x 5 -4 l 2 x 3 + 7 l 3 x 4 -3 l 4 x 5 1 2l x 3 -1 l 2 x 4 + 1 2l 3 x 5 ⎤ ⎥ ⎥ ⎦ . ( 15
)
Data required to build the TITOP model G f,pol (s) defined in Eq. ( 12) are obtained using following decompositions:

M rr M T rf M rf M ff = ρSl 55440 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 55440 27720l 462l 2 27720 -5544l 462l 2 27720l 18480l 2 198l 3 19800l -3432l 2 264l 3 462l 2 198l 3 6l 4 181l 2 -52l 3 5l 4 27720 19800l 181l 2 21720 -3732l 281l 2 -5544l -3432l 2 -52l 3 -3732l 832l 2 -69l 3 462l 2 264l 3 5l 4 281l 2 -69l 3 6l 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 16
)
K = EI z 70l 3 ⎡ ⎢ ⎢ ⎣ 6l 4 -30l 2 8l 3 l 4 -30l 2 1200 -600l 30l 2 8l 3 -600l 384l 2 -22l 3 l 4 30l 2 -22l 3 6l 4 ⎤ ⎥ ⎥ ⎦ , τ T Φ T (l) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 l 0 0 0 1 0 0 1 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 17
)
The clamped-free and free-clamped models obtained from the TITOP models G f,pol (s) and G f,n (s) (for various numbers n of assumed modes) are now evaluated in Table 1 regarding the accuracy of the first two natural frequencies. Of course, the model G f,n (s) provides the exact clamped-free frequencies for any values of n since the clamped-free modes are exactly taken into account in this approach. But its inverse, which must represents the freeclamped model, is not accurate at all, even for high values of n (for n = 16, the order of the model G f,16 (s) is 32!!). In comparison, the model G f,pol (s) is an eighth-order model with only four flexible modes, which is quite accurate regarding the first two natural frequencies from the clamped-free to free-clamped conditions. In the following section, the TITOP model G f,pol (s) is used to model an Euler-Bernoulli beam under the various well-known boundary conditions: free, clamped, support, or sliding on each of the two tips of the beam. In each case, the first two natural frequencies are compared with the reference values provided by the beam theory.

Table 1 Comparison of the first two frequencies of TITOP models and their inverses. The direct model corresponds to the clamped-free model of an Euler-Bernoulli beam. The inverse model corresponds to the free-clamped model. Frequencies are expressed using normalized units

( EI z ρSl 4 ) i ω i,ref G f,pol (s) G -1 f,pol (s) G f,n (s) G -1 f,4 (s) G -1 f,8 (s) G -1 f,

Other boundary conditions

The various classical (clamped, free, pinned, sliding) boundary conditions are taken into account by setting to zero the corresponding inputs of G f,pol (s), or of its inversion around one or more channels. The procedure for such a channel inversion is described in Appendix 1 under state-space formalism and can be directly applied to the TITOP model G f,pol (s) defined in Eq. ( 12). Note that the eight-component state vector x f of G f,pol (s) or of its various inverses has a physical meaning and is defined by

x f = q q ,
where more particularly (from ( 5) and ( 14))

x f (2) = q(2) = ṽ(l) = v C -v P -lα P , x f (3) = q(3) = ṽ (l) = α c -α P . ( 18
)
In the various cases presented in Appendix 2, the channel inversion to be used on the TITOP model G f,pol (s) in order to take into account the corresponding boundary conditions is detailed. Some boundary conditions can introduce also some internal state conditions, which can be used to reduce the order of the model. These conditions are also detailed in the various cases presented in Appendix 2. The first two natural frequencies ω i , i = 1, 2, are then compared with the reference values ω i,ref provided by beam theory [START_REF] Bishop | The Mechanics of Vibration[END_REF] under the same boundary conditions, and the relative errors ω i = 100

|ω i -ω ref i | ω ref i , i = 1, 2,
are computed and summarized in Tables 3 to 11.

Finally, we have to keep in mind that the model G f,pol (s) and its various inverses describe the relationship between forces and accelerations. These models provide the frequencies of the flexible modes for various boundary conditions. To take into account the rigid modes that can appear in some cases (free-free, pinned-free, sliding-free, sliding-sliding), we have to consider the positions v P , α P , v C , and α C as outputs of an augmented model. Of course, these positions can be obtained from accelerations using twice integrations, but it is recommended to use the internal state x f and Eqs. [START_REF] Perez | Flexible multibody system linear modeling for control using component modes synthesis and double-port approach[END_REF] to describe this augmented model by a minimal state-space representation. The block-diagram representation of this augmented model and the number of rigid modes are also detailed for the various cases presented in Appendix 2.

In every case, the TITOP model G f,pol (s) and its inverses provide a quite good approximation of the first two natural frequencies.

In conclusion, the polynomial approximation done in the TITOP model G f,pol (s) leads to a simple analytical model (order 8), which can be used during preliminary design phase.

It allows varying conditions on a very large range, from clamped to free conditions, for each translational or rotational d.o.f. The fifth-order polynomial shape functions chosen to build the TITOP model allow the first two natural frequencies to be representative enough for arbitrary boundary conditions. The only task needed to vary these boundary conditions consists of the inversion of input/output channels of the TITOP model or the feedback of some dynamic model between its input and output ports. That will be illustrated in Sect. 4.

Modal shapes of TITOP models

Let us consider the eighth-order TITOP model of the beam G f,pol (s) or its inverses G

-1 I f,pol (s).
For all these models, the state vector is still

x f (t) = [q T (t) qT (t)] T .
Let A, B, C, and D be the four state-space matrices of one of these models. Then the eigenvalue decomposition of A characterizes the free response of the model to initial conditions x f (0) (i.e., when the inputs are set to 0). The eigenvalues λ i of A are autoconjugate pairs associated with each natural frequency (λ i = ±jω i ) of the ith flexible mode (i = 1, . . . , 4). The associated eigenvector v i represents the magnitude of the ith mode on the components of the state vector x f :

x f (t) = q(t) q(t) = 4 i=1 g i v i e j (ω i t+ϕ i ) , ( 19
)
where the scalars g i and ϕ i depend on the initial conditions x f (0). From Eqs. ( 5) and ( 6), the deformation of the beam is

v(x, t) = 1 x φ T (x) ⎡ ⎣ v P (t) α P (t) q(t) ⎤ ⎦ ,
where φ(x) is defined by [START_REF] Pascal | Dynamics analysis of a system of hinge-connected flexible bodies[END_REF]. Modal shapes can be characterized by their contributions to the second time-derivative of the beam deflection v(x, t) (which is proportional to a factor -ω 2 i to the beam deflection v(x, t)). Indeed, let ψ i (x) be the modal shape of the ith mode; then v(x, t) can be decomposed as

v(x, t) = 4 i=1 g i ψ i (x)e j (ω i t+ϕ i ) and v(x, t) = - 4 i=1 g i ω 2 i ψ i (x)e j (ω i t+ϕ i ) = 1 x φ T (x) ⎡ ⎣ vP (t) αP (t) q(t) ⎤ ⎦ . ( 20
)
q(t) can be decomposed using [START_REF] Rong | Modified finite element transfer matrix method for eigenvalue problem of flexible structures[END_REF] and the state matrix A (see Appendix 1 for the submatrix notation):

q(t) = 4 i=1 g i A [5 : 8], : v i e j (ω i t+ϕ i ) = - 4 i=1 g i ω 2 i v i [1 : 4] e j (ω i t+ϕ i ) . ( 21
)
Fig. 9 Modal shapes of the first (top) and second (bottom) free-clamped modes: ψ i (x) (solid bold), reference φ i (x) (dashed), and error

φ i (x) -ψ i (x) (solid)
If the beam is clamped at P , then vP (t) = αP (t) = 0, and the identification of Eq. ( 20) with [START_REF] Schiehlen | Multibody system dynamics: roots and perspectives[END_REF] gives

ψ i (x) = -φ T (x)A [5 : 8], : v i /ω 2 i = φ T (x)v i [1 : 4] .
In the other cases, according to the channels (defined in the vector of indexes I) which are inverted in the model G -1 I f,pol (s), vP (t) and αP (t) can be outputs in the model, and their contribution to the modal shape ψ i (x) can be characterized using the matrix C of the model state-space representation. In a general way, we can define:

C v P = C(3, :) if 3 ∈ I, C v P = 0 1×8 otherwise, C α P = C(4, :) if 4 ∈ I, C α P = 0 1×8 otherwise.
Then,

ψ i (x) = -1 x φ T (x) ⎡ ⎣ C v P C α P A([5 : 8], :) ⎤ ⎦ v i /ω 2 i . ( 22
)
Thus, the modal shapes for arbitrary boundary conditions can be easily characterized from the eigenvectors/eigenvalues decomposition of the TITOP model or its inverses. As an example, the modal shapes ψ 1 (x) and ψ 2 (x) of the first two modes of the free-clamped beam modeled by G -1 [1:4] f,pol (s) are depicted in Fig. 9 and compared with the reference modal shapes φ 1 (x) and φ 2 (x) provided by Eq. ( 13) by changing x by lx. We can appreciate the accuracy of the first two modal shapes provided by the TITOP model.

Pure traction along (P 0 , x)-axis

Only one flexible mode is considered to represent the traction-compression along the (P 0 , x)-axis. Let u(x, t) be the axial deformation at any point of abscissa x (x ∈ [0; l]) due to axial forces -F x P and F x C applied to the beam at points P and C, respectively. Then, we assume that

u(x, t) = u P (t) + x l u C (t) -u P (t) = u P (t) + x l δ u (t) = 1 x l u P (t) δ u (t) ,
that is, the deformation at x is proportional to the total deformation δ u (t) of the beam. Then, the kinetic energy reads

T x = 1 2 l 0 ρS ∂u ∂t (x, t) 2 dx = 1 2 [ uP δu ]M x uP δu with M x = ρSl 1 1/2
1/2 1/3 , and the potential energy reads

V x = 1 2 l 0 ES ∂u ∂x (x, t) 2 dx = 1 2 [u p δ u ]K x u p δ u with K x = ES l 0 0 0 1 . The TITOP model G t (s) relative to the traction/compression dynamics along the (P 0 , x)- axis is ⎡ ⎢ ⎢ ⎣ δu δu üC F x P ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 0 1 0 0 -3E ρl 2 0 3 ρSl -3 2 -3E ρl 2 0 3 ρSl -1 2 3ES 2l 0 -1 2 -ρSl 4 ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣ δ u δu F x C üP ⎤ ⎥ ⎥ ⎥ ⎦ . ( 23
)

Two-port model G(s) of a beam in the planar case

From Eqs. ( 12) and ( 23), the 6 × 6 tenth-order TITOP model G(s) of the beam in the plane (P 0 , x, y) reads

G(s) = T G f (s) 0 4×2 0 2×4 G t (s) T T , ( 24
)
where the permutation matrix T is

T = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .
The block-diagram of G(s) is depicted in Fig. 10.

In [START_REF] Murali | Mechanical-attitude controller co-design of large flexible space structures[END_REF], the TITOP model of a uniform beam in detailed in the three-dimensional case. This parametric model is then used to model and design a boom linking a flexible deployable antenna to a flexible spacecraft. Each link is a uniform beam with mass density ρ = 2714 Kg/m 3 and modulus of elasticity E = 7.110 10 N/m 2 . The lumped masses of bearing assemblies in revolute joints between link 2 and link 3 and between link 3 and link 4 are equal: m l = 0.042 kg

Illustration: planar flexible four-bar mechanism

In this section, the two-input-two-output model approach presented in the previous section is applied to the well-known four-bar mechanism depicted in Fig. 11, taking into account the flexibility in the crank, in the coupler, and in the follower. This mechanism is a closedloop kinematic chain with one rigid degree-of-freedom and was largely studied in the past decades [START_REF] Kitis | Natural frequencies and mode shapes of flexible mechanisms by a transfer matrix method[END_REF][START_REF] Shigley | Theory of Machines and Mechanisms[END_REF][START_REF] Turcic | Dynamic analysis of elastic mechanism systems. Part I: applications[END_REF] and more recently [START_REF] Sitti | Piezoelectrically actuated four-bar mechanism with two flexible links for micromechanical flying insect thorax[END_REF] to design an efficient vibrating actuator. The objective is to compute the linear model G 4 bars (s, θ 1 ) between the torque T z P 1 applied on the crank at point P 1 and the angular acceleration of the crank αP 1 = δ θ1 for a given angular configuration θ 1 . First natural frequencies and modal shapes of the mechanism are then analyzed with respect to the crank angle θ 1 . The numerical values are taken from [START_REF] Kitis | Natural frequencies and mode shapes of flexible mechanisms by a transfer matrix method[END_REF] (translated in the International System of Units) for result comparison and summarized in Table 2. The corresponding natural frequencies are then called free frequencies (and denoted ω f (i)) since the crank is free to rotate around the (P 1 , z 1 ) axis. The modeling approach provides also the cantilevered frequencies ω c (i), thats is, when the crank is clamped on the ground, which are easily computed by considering the inverse model G -1 4 bars (s, θ 1 ). From data of Table 2, TITOP models of the crank G 1 (s), the coupler G 2 (s), and the follower G 3 (s) can be built from Eqs. ( 12), ( 16), ( 17), [START_REF] Shigley | Theory of Machines and Mechanisms[END_REF], and [START_REF] Simeon | On Lagrange multipliers in flexible multibody dynamics[END_REF]. The ground (i = 0 or 4) is assumed to be rigid. The planar wrenches and acceleration twists on the inputs and outputs of each model G i (s) are expressed in the body frame R i . Considering the linear model around a given angular configuration (small variations), all the frames R i are assumed to be inertial frames. For a given crank angle θ 1 , two kinematic constraints hold:

l 1 cos θ 1 + l 2 cos (θ 1 + θ 2 ) + l 3 cos (θ 1 + θ 2 + θ 3 ) = l 0 , l 1 sin θ 1 + l 2 sin (θ 1 + θ 2 ) + l 3 sin (θ 1 + θ 2 + θ 3 ) = 0.
They are solved in θ 2 and θ 3 using the equations given by Shigley and Uicker [START_REF] Shigley | Theory of Machines and Mechanisms[END_REF]. The values of θ i , i = 1, 2, 3, allow the (planar) direction cosine matrices between the various frames R i to be computed:

P i/i-1 = cos θ i -sin θ i sin θ i cos θ i .
The revolute joints between the links enforce geometrical constraints between inputs and outputs of models G i (s):

P i ≡ C i-1 ⇒ üP i vP i = P T i/i-1 üC i-1 vC i-1 .
Some dynamic boundary conditions also hold on the inputs of the three models using channel(s) inversion if required:

-for the crank (i = 1) and the coupler (i = 2), the inputs of the model i must take into account that: (i) F x C i and F y C i are constrained by the contact force applied by link i + 1 and the inertial force of the lumped mass:

F x C i F y C i = P i+1/i F x P i+1 F y P i+1 -m l üC i vC i , (25)
(ii) T z C i = 0 due to the revolute joint (no torque can be applied though the passive revolute joint), (iii) -T z P i is equal to the torque u applied by the driven system (crank case) or T z P i = 0 if the revolute joint is free (coupler case). Thus, T z P i must be on the input, and the models to be used for the crank and the coupler are G -1 6 1 (s) and G -1 6 2 (s), respectively, -for the follower (i = 3), the same condition on T z P 3 appears, but it is necessary to invert the first two channels of G 3 (s) in order to close the kinematic loop, and to impose that üC 3 and vC 3 are zero due to the revolute joint on ground. Then, the model to be used for the follower is G All these geometrical and dynamical boundary conditions can be easily expressed using the block diagram model depicted in Fig. 12. From this block diagram representation, the model G 4 bars (s, θ 1 ) between the torque u applied to the crank and the crank acceleration y = αP 1 can then be easily derived using MATLAB/SIMULINK and the function linmod. The modal analysis of this model provides the free pulsations ω f (i) of the mechanism. Note that the order of G 4 bars (s, θ 1 ) is 30: this model includes 15 flexible modes (five modes per link with four modes for bending in the plane (P i , y i , z i ) and one mode for traction-compression along the (P i , x i )-axis). But this model does not include the rigid mode of the mechanism. For control design purpose, this rigid mode can be included augmenting the model by a double integration:

α p 1 u (s, θ 1 ) = G 4 bars (s, θ 1 ) 1 s 2 .
To evaluate the accuracy of this model, the cantilevered pulsations ω c (i) are computed from the inverse model G -1 4 bars (s, θ 1 ) and compared with the results already published on the same mechanism [START_REF] Kitis | Natural frequencies and mode shapes of flexible mechanisms by a transfer matrix method[END_REF], [START_REF] Turcic | Dynamic analysis of elastic mechanism systems. Part I: applications[END_REF]. The first, second, and third cantilevered frequencies are plotted in Fig. 13 versus crank angle. The modal shapes of the mechanism computed from the accelerations (see Sect. 3.1.3) of the first three modes for two crank angle configurations (0 and 180 • ) are plotted in Figs. 14 and15. We can see that these results fit exactly the ones presented in [START_REF] Kitis | Natural frequencies and mode shapes of flexible mechanisms by a transfer matrix method[END_REF].

Conclusions

The objective of this paper is to present the interest of the two-input-two-output port approach to model any kind of flexible multibody systems. The novelty lies in the manipulation of transfer between accelerations and forces and to embed, in the same minimal state-space representation, the direct dynamic model of a body seen from one point and the inverse dynamic model seen from another point. All channels of this transfer are invertible allowing us to consider arbitrary boundary conditions at these points. By comparison with the Newton-Euler method, commonly used to model multibody system dynamics, the forward and backward recursions are directly solved by the feedback connections between the inputs and the outputs of the TITOP models of each body. Such a dynamic model is today restricted to the linear behavior of the multibody system and to loads lumped at some particular points Further works need to be performed to extend this approach to aeronautical engineering and, more particularly, to take into account distributed loads along each body of the system and nonlinear effects, which could appear when the rate of variation of the configuration Fig. [START_REF] Perez | Linear dynamic modeling of spacecraft with open-chain assembly of flexible bodies for ACS/structure co-design[END_REF] Inversion of the ith channel of G parameters is determinant. In this current state, this approach is fully applicable to space applications (large flexible orbital structures), but the proposed extensions are still required to apply it, for instance, to the helicopter rotor dynamic model.

Multichannel inversion

Let I be the vector (with q components) of indices corresponding to the channels to be inverted. The successive inversion of the q channels in G(s) is denoted

G -1 I (s) = G -1 I(1) -1 I(2) ••• -1 I(q) (s). ( 28
)
Appendix 2: TITOP model of the EULER-BERNOULLI beam under various boundary conditions Case 1: free-free beam (Fig. 17) The model G -1 [START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF][START_REF] De Luca | Closed-form dynamic model of planar multilink lightweight robots[END_REF] f,pol (s) is used to set to zero its inputs F Case 2: clamped-clamped beam (Fig. 18) The model G -1 [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF][START_REF] Bishop | The Mechanics of Vibration[END_REF] f,pol (s) is used to set to zero its inputs vC , αC , vP , and αP . There are no rigid modes in this case. The internal state constraints

ṽ(l) = x f (2) = 0, ṽ (l) = x f (3) = 0 ∀t ⇒ ẋf (2) = ẋf (3) = ẋf (6) = ẋf (7) = 0
can be used to truncate the state-space representation of G -1 [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF][START_REF] Bishop | The Mechanics of Vibration[END_REF] f,pol (s) to a fourth-order model. The frequencies of the two flexible modes are displayed in Table 4. (s) is used to set to zero its inputs vC , αC , vP , and T z P . There are no rigid modes in this case. The internal state constraints ṽ(l)l ṽ (l) = x f (2)lx f (3) = 0 ∀t ⇒ ẋf (2)l ẋf (3) = ẋf (6)l ẋf (7) = 0 can be used to reduce the state-space representation of G -1 [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF][START_REF] Bishop | The Mechanics of Vibration[END_REF][START_REF] De Luca | Closed-form dynamic model of planar multilink lightweight robots[END_REF] f,pol (s) to a sixth-order model. Case 5: pinned-pinned beam (Fig. 21) The model G -1 [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF][START_REF] De Luca | Closed-form dynamic model of planar multilink lightweight robots[END_REF] f,pol (s) is used to set to 0 its inputs: vC , T z c , vP , and T z P . There are no rigid modes in this case. Case 7: sliding-clamped beam (see Fig. 23) The model G -1 [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF][START_REF] Bishop | The Mechanics of Vibration[END_REF][START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF] f,pol (s) is used to set to zero its inputs vC , αC , F y P , and αP . There are no rigid modes in this case. The internal state constraints ṽ (l) = x f (3) = 0 ∀t ⇒ ẋf (3) = ẋf (7) = 0 can be used to reduce the state-space representation of G -1 [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF][START_REF] Bishop | The Mechanics of Vibration[END_REF][START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF] f,pol (s) to a sixth-order model. Case 8: sliding-pinned beam (see Fig. 24) The model G -1 [START_REF] Alazard | Two-input two-output port model for mechanical systems[END_REF][START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF] f,pol (s) is used to set to zero its inputs vC , T z C , F y P , and αP . There are no rigid modes in this case.

Fig. 3

 3 Fig. 3 Block-diagram of the TITOP model G

Fig. 4 A

 4 Fig. 4 A flexible boom B linking a flexible antenna A to a flexible spacecraft S through a revolute joint (left) and its block diagram model (right)

Fig. 5 AFig. 6 A

 56 Fig. 5 A flexible closed-loop structure (left) and its block diagram model (right)

Fig. 7

 7 Fig. 7 Parameterization of the beam deflection at a given time t

0

  xφ(x) dx] is the n × 2 rigid-flexible coupling mass matrix, -M ff = ρS l 0 φ(x)φ(x) T dx is the n × n flexible mass matrix.

10 ) with τ = 1 l 0 1 ,

 101 the rigid kinematic model between point C and P of the beam in the planar case, and Φ(l) = [φ(l) φ (l)] T . The relation between accelerations at points C and P , under the low-speed motion assumption reads vC αC = τ vP αP + Φ(l) q.

Fig. 8

 8 Fig. 8 TITOP model G f (s) block diagram of a flexible beam in bending

  n, are the n first solutions of the characteristic equation cos βl cosh βl + 1 = 0. The natural frequencies are given by ω i,ref = β 2 i EIz ρS and are the clamped-free reference values to evaluate approximate solutions. Such a decomposition basis is orthonormal and leads to M ff = I n and K = diag([ω 2 i,ref ]). The values of M rf and Φ(l) are not detailed for brevity.

Fig. 10 Fig. 11

 1011 Fig. 10 TITOP model G(s) block diagram

3 (

 3 s) and the first two outputs are the forces F x C 3 and F y C 3 applied by the follower on the ground.

Fig. 12

 12 Fig. 12 Block-diagram of model G 4 bars (s, θ 1 ) between the crank torque u and the crank acceleration y = αP 1

Fig. 13 1 Fig. 14 Fig. 15

 1311415 Fig. 13 ω c (i), i = 1, 2, 3, versus crank angle θ 1

  T z P . The augmented model with the two rigid modes is depicted in Fig.17

  (right).

Fig. 17

 17 Fig.[START_REF] Perez | A flexible appendage model for use in integrated control/structure spacecraft design[END_REF] The free-free beam (left) and its block-diagram model (right)

Fig. 18

 18 Fig. 18 The clamped-clamped beam (left) and its block-diagram model (right)

Case 3 :Fig. 19 ) 4 f

 3194 pinned-free beam (The model G -1 ,pol (s) is used to set to zero its inputs F y C , T z C , vP , and T z P . The augmented model with the rigid mode is depicted in Fig. 19 (right).

Fig. 19

 19 Fig.[START_REF] Rong | Modified finite element transfer matrix method for eigenvalue problem of flexible structures[END_REF] The pinned-free beam (left) and its block-diagram model (right)

Fig. 20

 20 Fig.[START_REF] Rui | Transfer matrix method for linear multibody system[END_REF] The pinned-clamped beam (right) and its block-diagram model (left)

Fig. 21

 21 Fig.[START_REF] Schiehlen | Multibody system dynamics: roots and perspectives[END_REF] The pinned-pinned beam (left) and its block-diagram model (right)

Case 6 :

 6 sliding-free beam (Fig. 22) The model G -1 3 f,pol (s) is used to set to zero its inputs F y c , T z c , F y P , and αP . The augmented model with the rigid mode is depicted in Fig. 22 (right).

Fig. 22

 22 Fig.[START_REF] Shabana | Flexible multibody dynamics: review of past and recent developments[END_REF] The sliding-free beam (left) and its block-diagram model (right)

Fig. 23

 23 Fig.[START_REF] Shigley | Theory of Machines and Mechanisms[END_REF] The sliding-clamped beam (right) and its block-diagram model (left)

Fig. 24

 24 Fig.[START_REF] Simeon | On Lagrange multipliers in flexible multibody dynamics[END_REF] The sliding-pinned beam (right) and its block-diagram model (left)

Case 9 :Fig. 25 )

 925 sliding-sliding beam (The model G -1 [2,3] f,pol (s) is used to set to zero its inputs F y c , αc , F y P , and αP . The augmented model with the rigid mode is depicted in Fig. 25 (right).

Fig. 25

 25 Fig.[START_REF] Sitti | Piezoelectrically actuated four-bar mechanism with two flexible links for micromechanical flying insect thorax[END_REF] The sliding-sliding beam (left) and its block-diagram model (right)

Table 2

 2 Four-bar mechanism links properties

	i	0	1	2	3
	Name	ground	crank	coupler	follower
	Length l i (m)	0.254	0.108	0.2794	0.2705
	Cross section S i (m 2 )	-	1 .0774 × 10 -4	4.0645 × 10 -5	4.0645 × 10 -5
	Flexural rigidity (EI ) i (Nm 2 )	-	1 1 .472	0.616	0.616

Table 3

 3 The first two natural frequencies ω i of the free-free

	beam model G f,pol (s) and the -1 [3,4]	i	ω i,ref (	EI z ρSl 4 )	ω i (	EI z ρSl 4 )	ω i (%)
	relative errors with the reference	1	2 2 .373		22.564	0.85
	values ω i,ref	2	6 1 .673		63.537	3.0

Table 4

 4 The first two natural frequencies ω i of the clamped-clamped beam model

		i	ω i,ref (	EI z ρSl 4 )	ω i (	EI z ρSl 4 )	ω i (%)
	G -1 [1,2] f,pol (s) and the relative	1	2 2 .373		22.45		0.34
	errors with the reference values ω i,ref	2	6 1 .673		62.929	2.0

Table 5

 5 The first two natural frequencies ω i of the pinned-free beam model G

	-1 4 f,pol (s) and the	i	ω i,ref (	EI z ρSl 4 )	ω i (	EI z ρSl 4 )	ω i (%)
	relative errors with the reference values ω i,ref	1	1 5 .418		15.445	0.17

Table 6

 6 The first two natural frequencies ω i of the pinned-clamped beam model

			i	ω i,ref (	EI z ρSl 4 )	ω i (	EI z ρSl 4 )	ω i (%)
	G -1 [1,2,4] f,pol	(s) and the relative	1	1 5 .418		15.433	0.1
	errors with the reference values ω i,ref	2	4 9 .965		50.724	1.5

Table 7

 7 The first two natural frequencies ω i of the pinned-pinned beam model

		i	ω i,ref (	EI z ρSl 4 )	ω i (	EI z ρSl 4 )	ω i (%)
	G -1 [1,4] f,pol (s) and the relative	1	9 .8696		9.8725	0.03
	errors with the reference values ω i,ref	2	3 9 .478		39.646	0.4

Table 8

 8 The first two natural frequencies ω i of the sliding-free beam model G

	f,pol (s) and the -1 3	i	ω i,ref (	EI z ρSl 4 )	ω i (	EI z ρSl 4 )	ω i (%)
	relative errors with the reference values ω i,ref	1 2	5 .5932 3 0 .226		5.5934 30.734	0.003 1.7

Table 9

 9 The first two natural frequencies ω i of the sliding-clamped beam model

			i	ω i,ref (	EI z ρSl 4 )	ω i (	EI z ρSl 4 )	ω i (%)
	G -1 [1,2,3] f,pol	(s) and the relative	1	5 .5932		5.5933	0.002
	errors with the reference values ω i,ref	2	3 0 .226		30.561	1.1

Table 10

 10 The first natural frequencies ω i of the sliding-pinned beam model

		i	ω i,ref (	EI z ρSl 4 )	ω i (	EI z ρSl 4 )	ω i (%)
	G -1 [1,3] f,pol (s) and the relative	1	2 .4674		2.4674	0.0
	errors with the reference values ω i,ref	2	2 2 .207		22.272	0.3

Table 11

 11 The first two natural frequencies ω i of the sliding-sliding beam model

		i	ω i,ref (	EI z ρSl 4 )	ω i (	EI z ρSl 4 )	ω i (%)
	G -1 [2,3] f,pol (s) and the relative	1	9 .8696		9.8697	0.0008
	errors with the reference values ω i,ref	2	3 9 .478		40.988	3.8

Acknowledgements This research is supported, in part, by ONERA (The French Aerospace Lab), CNES (The French Space Agency), ESA (European Space Agency), Thales Alenia Space and Polytechnique Montréal. We would like to thank more particularly Christelle Cumer, Thomas Loquen, Christelle Pittet, Finn Ankersen, Luca Massotti, Catherine le Peuvédic, Chiara Toglia, David-Alexandre Saussié, and Hari Murali.

Appendix 1: Matrix and linear system operations

In this section, we present basic operations on matrices and linear dynamic systems and finally the inversion of one (or more) input-output channel(s) in a linear system.

Let us consider a square (the same number of inputs and outputs) linear system G(s) (s stands for the Laplace variable) with order n and m channels (i.e., m inputs and m outputs). A state-space realization of the system G(s) is a set of four matrices A n×n , B n×m , C m×n , and D m×m such that

Submatrices and subsystems MATLAB nomenclature is used to define the submatrix M(I, J) (resp. subsystem G I,J (s)) restricted to the rows (resp. outputs) and the columns (resp. inputs) ordered in vectors of indices I and J.

Single-channel inversion

The system corresponding to the inversion of the ith channel of the system G(s)

and can be characterized by the following state-space realization (Eq. ( 27)). Let J be the vector of indices form 1 to m without i:

In G -1 i (s), the ith inverted channel appears on the last channel; it is thus required to reorder the channels using the vector of indices K = [1, . . . , i -1, m, i + 1, . . . , m -1], and then

Let u and y be the input and output vectors of G. This inversion can be represented by the block diagram depicted in Fig. 16.


Antonia Estero-Botaro

Antonio García-Domínguez
email: antonio.garciadominguez@uca.esjuanjose.dominguez

Juan José Domínguez-Jiménez

Francisco Palomo-Lozano

Inmaculada Medina-Bulo
email: inmaculada.medina@uca.es

A framework for Genetic Test-Case Generation for WS-BPEL Compositions

Search-based testing generates test cases by encoding an adequacy criterion as the fitness function that drives a search-based optimization algorithm. Genetic algorithms have been successfully applied in search-based testing: while most of them use adequacy criteria based on the structure of the program, some try to maximize the mutation score of the test suite. This work presents a genetic algorithm for generating a test suite for mutation testing. The algorithm adopts several features from existing bacteriological algorithms, using single test cases as individuals and keeping generated individuals in a memory. The algorithm can optionally use automated seeding when producing the first population, by taking into account interesting constants in the source code. We have implemented this algorithm in a framework and we have applied it to a WS-BPEL composition, measuring to which extent the genetic algorithm improves the initial random test suite. We compare our genetic algorithm, with and without automated seeding, to random testing.

Introduction

Search-based testing [START_REF] Harman | A theoretical and empirical study of search-based testing: Local, global, and hybrid search[END_REF] consists of generating test data according to a certain adequacy criterion, by encoding it as the fitness function that drives a search-based optimization algorithm. Evolutionary testing is a field of searchbased testing that uses evolutionary algorithms to guide the search. The global searches performed by these algorithms are usually, but not always, implemented as genetic algorithms (GAs) [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF][START_REF] Holland | Adaptation in Natural and Artificial Systems[END_REF].

Using GAs for generating test data dates back to the work by Xanthakis et al. [START_REF] Xanthakis | Application of genetic algorithms to software testing[END_REF]. Since then, various alternative approaches for generating test data have been proposed, both using GAs and other evolutionary techniques. Mantere and Alander [START_REF] Mantere | Evolutionary software engineering, a review[END_REF] published a review of the works that applied evolutionary algorithms to software testing at the time.

Search-based testing has been mostly used for structural test data generation, with branch coverage as the most common adequacy criterion. The survey by McMinn [START_REF] Mcminn | Search-based software test data generation: A survey[END_REF] classified evolutionary algorithms for structural test data generation by the type of information used by the fitness function.

Mutation testing [START_REF] Offutt | Mutation Testing for the New Century, chap[END_REF] is a testing technique injecting simple faults in the program under test through the use of mutation operators. As a result, we obtain mutants: variants of the original program. The original program and its mutants are executed on a given test-suite. When the output of a mutant for a test-case does not agree with the output of the original program for the same test-case, the mutant has been killed by that test-case and, so, it is dead. This means that the test-case has served a purpose: detecting the fault that is present in the mutant. If the output is always the same as the original program for every test case in the test-suite, then the mutant remains alive. When this is always the case, regardless the input, the mutant is said to be equivalent to the original program.

Mutation testing uses the mutation score to measure the quality of the testsuite. The mutation score is the ratio of the number of mutants killed by the test-suite to the number of non-equivalent mutants. Obtaining test-suites with mutation scores in the 50%-70% range can be relatively easy, but increasing these scores to the 90%-100% can be very hard. For this reason, developing new techniques for generating test-suites is currently a very active field of study [START_REF] Baudry | From genetic to bacteriological algorithms for mutation-based testing[END_REF].

Several evolutionary approaches for generating test data for mutation testing already exist. Baudry et al. [START_REF] Baudry | From genetic to bacteriological algorithms for mutation-based testing[END_REF] compared genetic and bacteriological algorithms for generating test data for a mutation testing system. Likewise, May et al. [START_REF] May | Immune and evolutionary approaches to software mutation testing[END_REF] compared genetic and immunological algorithms for the same purpose. In both cases, the alternative worked better than the GA.

This work presents an evolutionary technique for test-case generation using a GA for a mutation testing system. Our GA adopts some features from the bacteriological algorithm proposed by Baudry et al. We built a framework to apply it to WS-BPEL (Web Service Business Process Execution Language [START_REF]OASIS: Web Services Business Process Execution Language 2[END_REF]) compositions. WS-BPEL was targeted because testing service-oriented software is not an easy task, and this language is an OASIS standard and an industrialstrength alternative in its field. WS-BPEL compositions may build new WS from other WS from all around the world, while using unconventional and advanced programming concepts as well. As the relevance and economic impact of service compositions grow, the need for efficient and effective testing techniques for this kind of software increases.

This work is structured as follows. Section 2 introduces the WS-BPEL language and the basic concepts behind GAs. Section 3 discusses related work. Section 4 describes the overall design of the GA for generating test-suites. Section 5 describes how the GA can be applied to WS-BPEL compositions and shows the results obtained by the GA generating tests for a particular WS-BPEL composition. Finally, we present our conclusions and future work in Section 6.

Background

A Brief Introduction to WS-BPEL

WS-BPEL 2.0 is an OASIS standard [START_REF]OASIS: Web Services Business Process Execution Language 2[END_REF] and an industrial-strength programming language for WS compositions. WS-BPEL is an XML-based language al- lows one to specify the behavior of a business process based on its interactions with other WS. The major building blocks of a WS-BPEL process are activities. There are two types: basic and structured activities. Basic activities only perform one purpose (receiving a message from a partner, sending a message to a partner, assigning to a variable, etc.). Structured activities define the business logic and may contain other activities. Activities may have both attributes and a set of containers associated to them. Of course, these containers can include elements with their own attributes too. We can illustrate this with the example skeleton in Figure 1.

WS-BPEL provides concurrency and synchronization between activities. An example is the flow activity, which launches a set of activities in parallel and allows to specify the synchronization conditions between them. In the aforementioned example we can see a flow activity that invokes three WS in parallel: checkFlight, checkHotel, and checkRentCar. Moreover, there is another WS, bookFlight, that will be invoked upon checkFlight completion. This synchronization between activities is achieved by establishing a link, so that the target activity of the link will be eventually executed only after the source activity of the link has been completed.

Genetic Algorithms

Genetic algorithms [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF][START_REF] Holland | Adaptation in Natural and Artificial Systems[END_REF] are probabilistic search techniques based on the theory of evolution and natural selection proposed by Charles Darwin, which Herbert Spencer summarized as "survival of the fittest".

GAs work with a population of solutions, known as individuals, and process them in parallel. Throughout successive generations, GAs perform a selection process to improve the population, so they are ideal for optimization purposes. In this sense, GAs favor the best individuals and generate new ones through the recombination and mutation of information from existing ones. The strengths of GAs are their flexibility, simplicity and ability for hybridization. Among their weaknesses are their stochastic and heuristic nature, and the difficulties in handling restrictions.

There is no a single type of GA, but rather several families that mostly differ in how individuals are encoded (binary, floating point, permutation, . . . ), and how the population is renewed in each generation. Generational GAs [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF] replace the entire population in each generation. However, steady-state or incremental GAs [START_REF] Syswerda | A study of reproduction in generational and steady-state genetic algorithms[END_REF] replace only a few (one or two) members of the population in each generation. Finally, in parallel fine-grained GAs [START_REF] Mühlenbein | Evolution in time and space -the parallel genetic algorithm[END_REF] the population is distributed into different nodes.

As each individual represents a solution to the problem to be solved, its fitness measures the quality of this solution. The average population fitness will be maximized along the different generations produced by the algorithm. The encoding scheme and the individual fitness used are highly dependent on the problem to solve, and they are the only link between the GA and the problem [START_REF] Mitchell | An introduction to genetic algorithms[END_REF].

GAs use two types of operators: selection and reproduction. Selection operators select individuals in a population for reproduction attending to the following rule: the higher the fitness of an individual, the higher the probability of being selected. Reproduction operators generate the new individuals in the population: crossover operators generate two new individuals or offspring, from two preselected individuals, their parents. The offspring inherit part of the information stored in both parents. On the other hand, mutation operators aim to alter the information stored in a given individual. The design of these operators heavily depends on the encoding scheme used.

Please, notice that the above mutation operators are related to the GA and are different from those for mutation testing.

Related Work

Evolutionary algorithms in general and GAs in particular have been widely used in various fields of software engineering, and especially so in software testing. A large part of the works on search-based testing used structural criteria to generate test cases.

Bottaci [START_REF] Bottaci | A genetic algorithm fitness function for mutation testing[END_REF] was the first to design a fitness function for a GA that generated test data for a mutation testing system. This fitness function is based on the three conditions listed by Offutt [START_REF] Offutt | Automatic test data generation[END_REF]: reachability, necessity and sufficiency. Being the first of its kind, it has been used in many works, even with other evolutionary algorithms. As an example, Ayari [START_REF] Ayari | Automatic mutation test input data generation via ant colony[END_REF] used it with ant colonies.

Baudry et al. [START_REF] Baudry | From genetic to bacteriological algorithms for mutation-based testing[END_REF] have compared GAs and bacteriological algorithms for generating test cases for a mutation system. The main issues of GAs were slow convergence and the need to apply the mutation genetic operator more often.

The main difference between the bacteriological algorithm and the GA is that it stores the best individuals and that it does not use the crossover operator. In addition, individuals now represent individual test cases instead of entire test suites. This alternative approach is shown to converge much faster than the previous one, producing higher mutation scores in less generations.

Our approach is based on a GA that has been modified to obtain some of the advantages of a bacteriological algorithm. The individuals of our GA are individual test cases and both crossovers and mutations are performed. We do not need to use a separate store for the best individuals, as we store all individuals in a memory to preserve the useful information that may reside in individuals with low fitness scores. Our fitness function differs from that in Baudry et al., as it takes into account both the number of mutants killed by the test-case and how many other test cases killed those mutants.

May et al. [START_REF] May | Immune and evolutionary approaches to software mutation testing[END_REF] presented two evolutionary approaches for obtaining test cases for a mutation testing system: one based on an immune approach, and a GA. The algorithm based on the immune approach iteratively evolved a population of antibodies (individual test cases), looking for those that killed at least one mutant not killed by any previous antibody. The selected antibodies were added to the internal set which would become the final output of the process.

Approach

We advocate an approach based on a GA whose goal is generating test cases killing the mutants generated from the program under test. Figure 2 illustrates the architecture underlying the test generation procedure. The test-case generator consists of two main components: the preprocessor and the GA.

Original program Preprocessor

Mutant generator

Mutants Genetic algorithm

Genetic test case generator

Test cases

Fig. 2. Architecture of the test-case generator

The preprocessor can produce all the information required by the GA to generate the test cases from the original program, providing the GA with an initial population of randomly generated test cases, completed with some test cases generated by automated seeding. The actual operations run by the preprocessor depend on the particular programming language in which the program under test is written. Please see Section 5.1 for the specific details.

The test cases produced in each generation are executed against all the mutants. The fitness of each test-case is computed by comparing the behavior of each mutant and the original program on it. This guides the search of the GA. As the GA completes its execution, a final test-suite is obtained.

The goal of a test-case generator is producing an optimal test-suite. We propose a GA with a design based on several features of bacteriological algorithms. The individuals of our GA are individual test cases. The GA will gradually improve the initial test-suite, generation by generation, using a memory similar to a hall of fame. The memory will store all the test cases produced so far. Therefore, the output of the GA will not be its last generation, but rather all the test cases in the memory. This memory plays the same role as the memory used in the bacteriological algorithm proposed by Baudry et al. [START_REF] Baudry | From genetic to bacteriological algorithms for mutation-based testing[END_REF].

Next, we discuss the particular choices that characterize this GA, paying particular attention to how they contribute to produce good test cases.

Encoding of individuals Each individual encodes one test-case. As test cases are highly dependent on the program under test, the structure of individuals should be flexible enough to adapt to any program. In order to solve this problem, each individual will contain an array of key-value pairs, where the key is the type of a program variable and the value is a particular literal of the corresponding type. Figure 3 illustrates how individuals are encoded. Fitness of individuals Individual fitness is a function of the number of mutants killed by the corresponding test-case. However, though tempting, just taking the number of killed mutants per test-case is not an appropriate fitness metric. Let us consider specialized test-cases able to kill just one mutant which can not be killed by any other test-case in the population. A metric based only in the number of kills would assign this important test-cases a low fitness.

For this reason, an additional variable has to be taken into account when assessing the fitness of individuals: the number of test-cases killing the same mutant. This is key in both distinguishing individuals with the same number of kills and being fair in the evaluation of specialized test-cases.

Let M be the set of mutants generated from the original program and T be the set of test-cases in a given population. The fitness function is computed from the execution matrix, E = [e ij ], where e ij is the result of comparing the execution of the original program versus the execution of mutant m j ∈ M on test-case t i ∈ T : e ij = 0 when no difference can be found, while e ij = 1 if m j was killed by t i . We define the fitness of individual t i as:

f (t i ) = |M | j=1 e ij |T | k=1 e kj (1)
Initial population We have implemented two ways of generating the initial population at the beginning of the execution of the GA:

1. The initial population corresponds to valid random data. 2. The initial population corresponds to valid random data with additional test-cases where some components are replaced by constants of matching types found in the source code. This approach is similar to those proposed in [START_REF] Alshahwan | Automated web application testing using search based software engineering[END_REF], [START_REF] Alshraideh | Search-based software test data generation for string data using program-specific search operators[END_REF], and [START_REF] Fraser | The seed is strong: Seeding strategies in search-based software testing[END_REF]. We have named this approach automated seeding.

Generations The test-case generator is based on a generational GA in which the population of the next generation consists of offspring produced from pairs of individuals of the population of the previous generation by the genetic operations of crossover and mutation.

The selection scheme determines the way that individuals are chosen for haploid reproduction, and eventually mutation, depending on their fitness. We select individuals with the roulette wheel method designed by Goldberg [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF]. Therefore, we make selection probability proportional to fitness.

Algorithm 1 illustrates how crossover and mutation operators can be applied to generate a pair of offspring from the current population. We assume that p c is the crossover probability, that p m is the probability of mutation, and that random-uniform(a, b) produces a random number uniformly distributed in the real interval [a, b). This procedure is iterated until a new population of identical size is reached. New populations will be generated until a given termination condition is met.

In a whole population of size n, crossover contributes a total of np c individuals. As mutation does not alter the number of individuals at all, n(1p c ) selected individuals remain unchanged to maintain the population size constant.

Crossover operator

The crossover genetic operator exchanges the components of parents to produce two offspring. We use one-point crossover, in which the point is randomly chosen using a uniform distribution. Figure 4 shows how crossover is performed on a pair of individuals.

Selection with Goldberg's roulette-wheel method. parents ← select-parents(population)

Crossover with probability pc. if random-uniform(0, 1) < pc then offspring ← do-crossover (parents) else offspring ← parents end if

Mutation with probability pm. if random-uniform(0, 1) < pm then offspring ← do-mutation(offspring) end if return offspring Algorithm 1: Offspring generation. int 8000 string "true" string "high" int 1500 string "true" string "low" int 8000 string "true" string "low" int 1500 string "true" string "high" Fig. 4. Application of the crossover operator Mutation operator The mutation genetic operator changes the value field of one component in an individual. The actual change will depend on its data type. For instance, a float value v is mutated into v using the following formula:

v = v + s • r p m ( 2
)
where s is randomly chosen between -1 or 1 and r is a parameter which modulates the change of v. For instance, if r = 1000, v may be incremented or decremented by up to 1000/p m . To summarize: v undergoes a perturbation that may be positive or negative, depending on the value of s. This perturbation is proportional to the inverse of p m : as mutations become less common, perturbations are larger.

Values of type string are mutated by replacing them with a new value.

When mutating list values, the operator randomly decides (with equal probabilities) whether to mutate its length or mutate an element picked at random. The length l of the list is mutated into a new length l as follows:

l = random-uniform(0, 1) • l p m (3)
Finally, tuple values are mutated by picking one of its elements at random, as their lengths are fixed.

Termination Termination conditions are checked on each generation after computing the fitness of the individuals. Our GA implements four termination conditions that can be combined or used in isolation: maximum number of generations, percentage of mutants killed, stagnation of maximum fitness, and stagnation of average fitness.

Application to WS-BPEL

In order to evaluate our approach we have implemented a framework for WS-BPEL compositions. By design, the framework encapsulates all the languagespecific (in this case, WS-BPEL-specific) details in a component known as the Preprocessor. Mutants are generated, executed, and evaluated using our mutation testing tool for WS-BPEL, MuBPEL [START_REF] Estero-Botaro | Quality metrics for mutation testing with applications to WS-BPEL compositions[END_REF], which incorporates the mutation operators defined by Estero et al. [START_REF] Estero-Botaro | Mutation operators for WS-BPEL 2.0[END_REF].

Preprocessor

A test-case for a WS-BPEL composition consists of the messages that need to be exchanged among the various partners in a business process: the client, the composition itself and the external services invoked by the composition. Generating test-cases automatically requires knowing the structure of the messages that constitute these test-cases.

In order to obtain the test-case file used by the GA as the initial population, an Analyzer produces a message catalog from the WSDL documents describing the public interfaces of the WS-BPEL composition. The message catalog contains a set of templates that can generate the required messages. Each template declares the variables used in it and their types. From this message catalog, a specification of the test data format is produced. This specification is used by the random test generator in the framework to produce a test data file.

We will use the LoanApproval composition [START_REF]OASIS: Web Services Business Process Execution Language 2[END_REF] to illustrate these steps. This composition simulates a loan-approval service in which customers request a loan of a certain amount. Whether the loan is approved or denied depends on the requested amount and the risk level associated to the customer. Two external WS are used: the assessor service and the approval service. When the requested amount is modest (less than or equal to 10000 monetary units) the assessor WS is invoked to assess the risk presented by the customer. If the risk is low, the loan is immediately approved. In any other case, that is, large loans or high-risk customers, the decision is delegated to the approval WS.

The LoanApproval composition uses three WSDL files: one for each of the external WS and another for the composition itself. From these files, the Analyzer obtains the message catalogs describing the variables from which the messages exchanged between the composition and its partners can be generated. The LoanApproval composition uses a single input variable named req_amount of type int. The approval service is invoked by the WS-BPEL composition, and its output is controlled by the Boolean variable ap_reply. The assessor service is also invoked by the composition, and its output variable as_reply accepts a string value between "high" or "low" which just represents the estimated risk.

Figure 5 shows the test data format specification that was extracted from the message catalog and that will be used to generate the test-cases. It is written in a domain-specific language used by the random test generator in the framework, which can produce files such as the one shown in Figure 6. This file shows seven test cases. Variable req_amount is set to 54907 in the first test case, 103324 in the second test-case and so on. Variables ap_reply and as_reply contain the replies to be sent from the mockups of the approver and the assessor services, respectively, and can take the values "true" or "false" and "high" or "low". The special value "silent" indicates that the mockup will not reply.

typedef int (min = 0, max = 200000) Quantity; typedef string (values = { "true", "false", "silent" }) ApReply; typedef string (values = { "low", "high", "silent" }) AsReply; Quantity req_amount; ApReply ap_reply; AsReply as_reply; = ["false", "silent", "false", "false", "true", "true", "silent"]) #set($as_reply = ["high", "silent", "low", "high", "low", "low", "high"])

Fig. 6. Test data file for the LoanApproval composition

Figure 7 shows a potential individual for the LoanApproval composition. We can see that a test-case for this composition consists of one component of type int and two string components that correspond to the variables req_amount, ap_reply and as_reply respectively. These were shown in Figure 6 as well.

int 8000 string "true" string "high" Fig. 7. Individual representing a test-case for the LoanApproval composition

A Case Study: The Loan Approval Composition

In this section, we will conduct several experiments to provide answers to the following research questions: RQ1. How much does the GA improve the quality of the initial test suite? The initial population of the GA is randomly generated. We use three quality metrics to compare the initial population with the test-suite produced by the GA: mutation score, sentence coverage and condition/decision coverage. RQ2. Are the improvements just a consequence of the greater size of the final testsuite? The GA is a test-case generator. As such, the test suite generated can be much greater than the initial test-suite. Thus, it could be more powerful just because of its size. We compare the test suites produced by the GA with randomly generated test-suites of the same size, using the same metrics employed in RQ1.

We will use the LoanApproval composition from Section 5.1 to answer the proposed research questions. This composition has 154 LOC, it generates 90 mutants, 2 mutants are invalid, and 9 mutants are equivalent. 1 We can observe that BPEL compositions tend to be smaller than traditional programs in the sense that they define the logic of the composition of the external WS or partners, while the bulk of the code is in the WS themselves. Therefore, the number of mutants obtained is comparatively much lower than in traditional programming languages.

Next we will describe the experimental procedure used to answer the research questions and summarize the results obtained.

Since the initial population of the GA is randomly generated, results could largely vary from one execution to another. In order to alleviate this, we decided to select as the initial population a good representative of possible initial populations through these steps:

1. Population sizes were defined. 2. Thirty random test-suites were generated, with sizes matching each population size. 3. Test-suites were run against the original composition and their mutants to measure their mutation scores, sentence coverages, and condition/decision coverages. 2 4. Medians for the quality metrics were computed. 5. A test-suite with median mutation score was selected.

The result of this process is a typical test-suite, which has been generated at random, but it is unbiased. Thus, it is a good candidate for initial population in our experiments. 1 Equivalent mutants have been manually detected by inspecting the source code of surviving mutants. 2 Sentence and condition/decision coverages have been computed as the percentage of killed mutants from those generated by a sentence coverage mutation operator and a condition/decision mutation operator, respectively.

Regarding population sizes, it is well known that a GA is likely to find better solutions when the population size is greater. But bigger population sizes demand more computational resources to find those solutions. However, Krishnakumar [START_REF] Krishnakumar | Microgenetic algorithms for stationary and nonstationary function optimization[END_REF] found that a population size of 10 individuals can produce similar convergence rates with a tailored GA.

Since the LoanApproval composition produces 90 mutants, a test suite might need up to 90 test-cases to kill them all, though much fewer test cases will usually suffice. For this reason, we defined population sizes in terms of percentages of the total number of mutants. Percentages of 15%, 20%, and 25% were selected as candidate values to produce small populations still having more than 10 individuals.

As to the number of executions, we decided to use 30 executions for each population size. Therefore, 30 different seeds were employed to generate random data for 30 different executions.

Table 1 shows the values of the three metrics for the selected test-suite (the random initial test-suite) and for the final test-suite produced by the GA, with and without automated seeding. These are the medians of all the values produced by the 30 executions for each population size. The configuration parameters for the GA are shown in Table 2. We set p m to 10%, as it is also suggested by Baudry et al. [START_REF] Baudry | From genetic to bacteriological algorithms for mutation-based testing[END_REF]. In agreement with Baudry's results, it was easy to obtain a mutation score within the 50%-70% range: in this particular case, the mutation score of the initial test-suite was 72%. Even without automated seeding, the three quality metrics have been considerably improved by the GA: mutation score increases from 0.72 to 0.96, and both coverage metrics increases from 75% to 100%. With automated seeding, mutation score further increases to 0.99. Since all the population sizes obtained the same medians in their quality metrics, we compared the results from each of their 30 executions without automated seeding and with automated seeding. The results are shown in Tables 3 and4, respectively. Without automated seeding, the size that produces the most stable results is 25%. With automated seeding, the size that produces the best results is 20%, reaching 100% mutation score in two executions. As the GA generates test-suites much bigger than the initial test-suites, the next step is comparing random generation of bigger test-suites with the results produced by our GA, with and without automated seeding. For a 20% population size, the GA executed a median of 156 distinct test-cases without automated seeding, and a median of 153 different test-cases with automated seeding.

Table 5 compares the median values of metrics corresponding to 30 random test-suites, with 156 test-cases each, with those produced by both variants of our GA. Results indicate that the metrics for random test-suites were as good as those for the GA without automated seeding. This can be explained by the fact that small compositions usually have few execution paths and, thus, random test-suites of considerable size are likely to produce high sentence and condition/decision coverages. Moreover, at the same time, those compositions tend to generate a small number of mutants. As a consequence, their mutants might be mostly killed by a random test-suite of even modest size and high mutation scores can be expected. This is the case, with the LoanApproval composition, which in fact has few mutants and execution paths. However, the GA with automated seeding obtained better mutation score (0.99) than the random test-suite (0.96). From these results, we can now answer the two research questions that were posed at the beginning of this section for the LoanApproval composition. RQ1. How much does the GA improve the quality of the initial test suite? The testsuite generated by our GA improves the initial test-suite with respect to all the quality metrics measured (mutation score, sentence coverage, and condition/decision coverage). Without automated seeding, the GA increases the mutation score from 0.72 to 0.96. Sentence coverage and condition/decision coverage increase from 75% to 100%. With automated seeding, the mutation score further increases to 0.99. RQ2. Are the improvements just a consequence of the greater size of the final testsuite? Once we extend the initial test-suite with sufficient random test-cases to match the size of the final test-suite, the GA shows no difference when automated seeding is disabled. However, with automated seeding, the GA obtains a higher mutation score: 0.99 instead of 0.96.

Conclusions and Future Work

We have presented a GA for generating test-cases for a mutation testing system. The initial population can be generated in two ways. The first way simply generates all individuals at random. The second way mixes in test-cases which are derived from the constants present in the source code of the program. Our GA also adopts several features of the bacteriological algorithms proposed for test-case generation by Baudry et al. [START_REF] Baudry | From genetic to bacteriological algorithms for mutation-based testing[END_REF].

We have implemented both approaches in a framework which generates testcases for WS-BPEL compositions. The GA was applied to a standard WS-BPEL composition. We compared the quality of the initial random test-suite against the test-suites produced by the GA. Test-suite quality was measured using mutation score, sentence coverage and condition/decision coverage. It is shown that, for this composition, the median of the quality metric values across 30 executions improves, both with and without automated seeding.

The quality of the test-suites generated by the GA has been also compared to the quality of random test-suites of the same size. In this case, random testing produces the same results as using the GA without automated seeding. However, the GA with automated seeding produces slightly better results. This has been traced to the fact that the composition under study is small and have few execution paths. Therefore, random test-suites of considerable size are likely to produce high sentence and condition/decision coverages. At the same time, its number of mutants is modest and they can be killed by a random test suite. Consequently, random test-suites enjoy higher mutation scores than usual and the margin of improvement for the GA, which starts from an initial population of random test-cases, gets drastically reduced.

These preliminary results are promising, but we acknowledge that one of the limitations of the present study is that it has only been applied to a single WS-BPEL composition so far. Future work will be devoted to validate our results with a number of increasingly complex WS-BPEL compositions. The main hindrance here is the absence of a standard set of WS-BPEL compositions suitable for functional testing.

Another limitation of the current study is that we have mainly concentrated in assessing the impact of the population size. Other configuration parameters, for which reasonable values were used, might be adjusted and hopefully improve our results. Further research and extensive experimentation could lead to better combinations of parameters for generating test-suites.

Fig. 1 .

 1 Fig. 1. WS-BPEL 2.0 activity sample

Type 1 Fig. 3 .

 13 Fig. 3. Encoding of individuals

Fig. 5 .

 5 Fig. 5. Specification of the data used to build messages for the LoanApproval composition

Table 1 .

 1 Quality metrics for the initial test-suite (random) and the test-suites generated by the GA. In columns: MS are mutation scores, SC are sentence coverages, and CDC are condition/decision coverages. Coverages are measured as percentages.

			Random			GA		GA (aut. seeding)
	Pop. size	MS	SC	CDC	MS	SC	CDC	MS	SC	CDC
	15%	0.72	75.0	75.0	0.96	100.0 100.0	0.99	100.0 100.0
	20%	0.72	75.0	75.0	0.96	100.0 100.0	0.99	100.0 100.0
	25%	0.72	75.0	75.0	0.96	100.0 100.0	0.99	100.0 100.0

Table 2 .

 2 Configuration parameters of the GA

	Population	pc	pm	r

Table 3 .

 3 Quality metrics obtained by the GA without automated seeding for each population size, grouped by execution results

	Population Seeds	MS	SC (%) CDC (%)
		2	0.72	75.0	75.0
	15%	1 1	0.90 0.95	93.8 100.0	100.0 100.0
		26	0.96	100.0	100.0
	20%	2 28	0.95 0.96	100.0 100.0	100.0 100.0
	25%	30	0.96	100.0	100.0

Table 4 .

 4 Quality metrics obtained by the GA with automated seeding for each population size, grouped by execution results

	Population Seeds	MS	SC (%) CDC (%)
	15%	7 23	0.97 0.99	100.0 100.0	100.0 100.0
		3	0.97	100.0	100.0
	20%	2	1.00	100.0	100.0
		25	0.99	100.0	100.0
	25%	3 27	0.97 0.99	100.0 100.0	100.0 100.0

Table 5 .

 5 Results produced by random test generation and the two versions of our GA

		Random GA GA (autom. seeding)
	Size test-suite	156	156	153
	MS	0.96	0.96	0.99
	SC (%)	100.0 100.0	100.0
	CDC (%)	100.0 100.0	100.0

Acknowledgments

This work was funded by the Spanish Ministry of Science and Innovation under the National Program for Research, Development and Innovation, project MoDSOA (TIN2011-27242).