Alexandr Klimchik

Stéphane Caro

Anatol Pashkevich

Anatol Pashkevich Practical

Practical identifiability of the manipulator link stiffness parameters

Keywords: Robotic machining, stiffness modeling, elastostatic identification, parameter identifiability

The paper addresses a problem of the manipulator stiffness modeling, which is extremely important for the precise manufacturing of contemporary aeronautic materials where the machining force causes significant compliance errors in the robot end-effector position. The main contributions are in the area of the elastostatic parameters identification. Particular attention is paid to the practical identifiability of the model parameters, which completely differs from the theoretical one that relies on the rank of the observation matrix only, without taking into account essential differences in the model parameter magnitudes and the measurement noise impact. This problem is relatively new in robotics and essentially differs from that arising in geometrical calibration. To solve the problem, several physical and statistical model reduction methods are proposed. They are based on the stiffness matrix sparseness taking into account the physical properties of the manipulator elements and also on the heuristic selection of the practically non-identifiable parameters that employs numerical analyses of the parameter estimates. The advantages of the developed approach are illustrated by an application example that deals with the stiffness modeling of an industrial robot used in aerospace industry.

INTRODUCTION

Application area of industrial robots is continuously increasing, they become more and more popular in many technological processes, including precise high-speed machining. For this process, the robot is subject to essential external loadings caused by the machining force that may lead to non-negligible end-effector deflections [START_REF] Dépincé | Active integration of tool deflection effects in end milling. Part 2. Compensation of tool deflection[END_REF]. This feature becomes extremely important in the aerospace industry, where the accuracy requirements are very high but the materials are hard to process. In this case, the manipulator stiffness modeling and corresponding error compensation technique are the key issues [START_REF] Kövecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF][START_REF] Eastwood | A gravitational deflection compensation strategy for HPKMs[END_REF][START_REF] Klimchik | Compliance error compensation technique for parallel robots composed of non-perfect serial chains[END_REF][START_REF] Chen | Spatial statistical analysis and compensation of machining errors for complex surfaces[END_REF].

In the manipulator stiffness modeling, there are three main approaches that can be summarized as follows: the Finite Element Analysis (FEA), the Matrix Structural Analysis (MSA), and the Virtual Joint Method (VJM). The most accurate of them is the FEA-based technique [START_REF] Piras | Dynamic finiteelement analysis of a planar high-speed, high-precision parallel manipulator with flexible links[END_REF], which allows presenting manipulator components with their true shape and dimension. However, this method is usually applied at the final design stage because of the high computational expenses [START_REF] Long | Optimal structural design of a planar parallel platform for machining[END_REF]. The MSA method [START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF] incorporates the main ideas of the FEA, but operates with rather large elements -3D flexible beams. This obviously leads to the reduction of the computational efforts, but does not eliminate the disadvantages of FEA. And finally, the VJM method [START_REF] Chen | Conservative Congruence Transformation for Joint and Cartesian Stiffness Matrices of Robotic Hands and Fingers[END_REF][START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF][START_REF] Quennouelle | Stiffness Matrix of Compliant Parallel Mechanisms[END_REF][START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF], is based on the extension of the traditional rigid model by adding the virtual joints (localized springs), which describe the elastic deformations of the links, joints and actuators. This technique provides reasonable trade-off between the model accuracy and computational complexity, which will be further used in this paper. However, evaluation or identification of the stiffness model parameters is not a trivial problem.

The main difficulty here is that the straightforward application of the VJM method yields very high number of parameters that differ in their magnitudes and in their impact on the model precision. Moreover, direct application of this technique may produce redundant models that are not suitable for calibration. In particular, attempts to solve the identification problem for the whole set of the elastostatic parameters (258 for 6 d.o.f. manipulator) leads to the fail of the numerical routines that is caused by singularity of the relevant observation matrix. It is worth mentioning that similar problem is also known in geometric calibration where the concept of completeirreducible-continues model has been introduced and relevant algebraic tools for the model reduction have been developed [START_REF] Everett | Kinematic modelling for robot calibration[END_REF][START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF][START_REF] Pashkevich | Computer-aided generation of complete irreducible models for robotic manipulators[END_REF]. However, in elastostatic calibration there is an additional difficulty caused by high number of relatively small parameters for which the measurement noise impact is very essential. As follows from our experience, the identification results may violate fundamental physical properties of the stiffness matrices, such as positive-definiteness and symmetry, and are not acceptable for the compliance error compensation (more details are given in Section 3.2 presenting a motivation example). For this reason, this paper introduces a new notion of practical identifiability and proposes corresponding model reduction methods that allow obtaining reliable results in real industrial environment.

To address the above mentioned problem, the remainder of the paper is organized as follows. Section 2 presents the stiffness modeling background. Section 3 describes the calibration procedure and contains the motivation example allowing us to define the research problems. In Section 4, the developed model reduction methods are presented. Section 5 contains an application example and illustrates advantages of the proposed technique. And finally, Section 6 summarizes the main contributions of the paper.

STIFFNESS MODELING BACKGROUND

Let us consider elastostatic model of a general serial manipulator, which consists of a fixed "Base", a serial chain of flexible "Links", a number of flexible actuated joints "Ac" and an "End-effector" (Fig. 1). It is assumed that all links are separated by either rotational or translational joints. Such architecture can be found most of industrial serial robots.

W W

Virtual springs

Figure 1 Serial manipulator and its VJM model

In order to evaluate the stiffness of the considered manipulator, let us apply the virtual joint method (VJM), which is based on the lump modeling approach [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. According to this approach, the original rigid model should be extended by adding virtual joints (localized springs), which describe elastic deformations of the links. Besides, virtual springs are included in the actuated joints, in order to take into account the stiffness of the control loop. Under these assumptions, the kinematic chain can be described by the following serial structure:

(a) a rigid link between the manipulator base and the first actuated joint described by the constant homogenous transformation matrix Base T ;

(b) several flexible actuated joints described by the homogeneous matrix function Joint Ac ( )

i i i q   T
, which depends on the actuated joint variable i q and the virtual joint variable

i i i i i i i      
 corresponding to the translation/rotation deflections in/around the axis x, y, z;

(e) a rigid link from the last joint to the end-effector, described by the constant homogenous matrix transformation Tool T .

In the frame of these notations, the final expression defining the end-effector location subject to variations of all joint coordinates may be presented as the product of the following homogenous matrices and matrix functions

  Base Joint Ac Link VJM Link Tool 1 () n i i i i i i q                T T T T T T ( 1
)
where n is the number of links/joints, and the components

Base Joint Link VJM Tool , , (.), , (.), i i
T T T T T may be factorized with respect to the terms including the joint variables (in order to simplify computing of the derivatives). For further convenience, after extraction from the homogeneous matrix T rotation and translation components [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods, and Algorithms[END_REF], the kinematic model can be rewritten in more conventional form ( , )  t g q

(2)

where (.) g denotes relevant vector function, the vector ( , , ..., ) q q q  q aggregates all actuated coordinates, the vector T 2 n θ 1 ( , , ..., )

  

 collects all virtual joint coordinates, and θ n is the number of the virtual joints. It should be noted that here the values of coordinates q are completely defined by the robot controller, while the values of the virtual joint coordinates θ θ θ

; 0        J t W J K (4)
describing elastostatic properties of the considered manipulator. In these equations, the end-effector displacement Δt is treated as the model input and the external wrench W is the model output, which corresponds to the representation of the manipulator stiffness matrix in the following form

C •Δ  W K t (5)
where C K is the desired Cartesian stiffness matrix of the considered manipulator for given robot configuration q . To find this matrix, equations (4) may be presented in the matrix form

θ T θ θ                 0J Wt 0 JK (6)
and solved for W . This transformation yields the following force-deflection relation

1 T θ θ θ Δ      J K J W
t that allows us to express the manipulator Cartesian stiffness matrix as

  1 1T C θ θ θ      K J K J (7)
This expression allows us to compute the Cartesian stiffness matrix assuming that the matrix

1 2 θ θ θ ( , , ...) diag  K K K
, defining elastostatic properties of the manipulator links/joins is given. However, in practice, the matrices  

(i)
θ , 1, 2,... i  K are unknown and should be identified from relevant experiments. However, there are a number of numerical problems that may arise here that are in the focus of the remaining parts of the paper.

PROBLEM OF ELASTOSTATIC PARAMETERS IDENTIFICATION

Methodology of elastostatic identification

To estimate the desired matrices describing elasticity of the manipulator components (i.e., compliances of the virtual springs presented in Fig. 1), the elastostatic model ( 5) should be rewritten as

  (i) (i) (i)T θ θ θ 1 • n i    J k J t W ( 8
)
where t is the vector of the end-effector displacements under the loading W , the matrices (i)

(i) 1 θ θ ()   kK
denote the link/joint compliances that should be identified via calibration, and the matrices (i) θ J are corresponding sub-Jacobians obtained by the fractioning of the aggregated Jacobian [START_REF] Dépincé | Active integration of tool deflection effects in end milling. Part 2. Compensation of tool deflection[END_REF] (2)

θ θ θ [ , ,...] T T T  J J J
. For the identification purposes, this expression should be transformed into more convenient form, where all desired parameters (elements of the matrices

(i) θ , 1, 2,... i  k
) are collected in a single vector [START_REF] Dépincé | Active integration of tool deflection effects in end milling. Part 2. Compensation of tool deflection[END_REF] (1) ( ) θ11 θ12 θ66

( , ,... )

n k k k  π . It yields the following linear equation )• (,    A q W t π (9)
where (.)

 A is so-called observation matrix that defines the mapping between the unknown compliances π and the endeffector displacements t under the loading W for the manipulator configuration q .

Taking into account that the calibration experiments are carried out for several manipulator configurations defined by the actuated joint coordinates , 1, j jm  q , the system of basic equations for the identification can be presented in the following form

( ; 1, • ,) jj j j j jm      t πε A q W ( 10
)
where j ε denotes the vector of measurement errors. Further, using these notations and assigning proper weights for each equation, the identification can be reduced to the following optimization problem

1 ( ) ( ) min m TT j j j j j F          π A π t A t π (11)
where is the matrix of weighting coefficients that normalizes the measurement data. This minimization problem yields the following solution

    1 1 1 ˆ• m m T T T T j j j j j j            π A A A t . ( 12
)
If the measurement noise is Gaussian (as it is assumed in conventional calibration techniques), expression [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF] provides us with a unbiased estimates for which

  Ê  ππ .
Corresponding covariance matrix evaluating the dispersion of the parameter estimate π from one identification session to another can be computed as follows

    1 2 1 1 1 1 ĉov( ) m m T T T T T j j j j j j m TT jj j                π A A A Σ A A A (13)
where the matrix

2 ( • ) T E  Σ ε ε
describes the statistical properties of the measurement errors.

It can be proved [START_REF] Rice | Mathematical Statistics and Data Analysis[END_REF] that the best results in terms of the identification accuracy are achieved if 1   Σ . It leads to the following covariance matrix of the manipulator compliance parameters

  1 1 1 ĉov( ) m TT j j j       π A Σ Σ A (14)
Such assignment of the weighing coefficients also allows us to avoid the problem of different units in the objective function [START_REF] Quennouelle | Stiffness Matrix of Compliant Parallel Mechanisms[END_REF], which arises in straightforward application of the leassquare technique to the robot parameters identification if the measurement system provides both position and orientation data. It should be noted that this particularity is usually omitted in conventional robot calibration. Another way to improve the identification accuracy is related to the proper selection of manipulator measurement configurations   , 1, j jm  q that is also known as the calibration experiment planning [START_REF] Daney | Choosing measurement poses for robot calibration with the local convergence method and Tabu search[END_REF], which directly influences on the observation matrices ( , ) jj j  A q W and on the covariance matrix [START_REF] Everett | Kinematic modelling for robot calibration[END_REF].

It is clear that expression [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF] gives reliable estimates of the parameters π if and only if the matrix

1 1 m T T j j j     A Σ Σ A
is invertible. It leads to the problem of the parameter identifiability that have been studied by a number of authors for the problem of geometrical calibration [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF][START_REF] Pashkevich | Computer-aided generation of complete irreducible models for robotic manipulators[END_REF]. Relevant techniques are based on the information matrix rank analysis (via either SVD-or QR-decomposition). However, in real industrial practice where the measurement not is non-negligible, the identifiable parameters are not equivalent in terms of accuracy (both absolute and relative) and expression [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF] can give rather surprising results for some of them. This motivates revision of the above mentioned notion (parameter identifiability) and its extension taking into account the identification accuracy defined by the covariance matrix [START_REF] Everett | Kinematic modelling for robot calibration[END_REF]. In the following sub-sections, the notion of practical identifiability is introduced and a motivation example is presented, which illustrates potential problems that may arises in the manipulator elastostatic calibrations if conventional techniques are applied.

Motivation example

To illustrate the problems that may arise in identification of the manipulator elastostatic parameters, let us consider a numerical example that deals with a single link of the Orthoglide manipulator (Fig. 2). Its compliance matrix has been obtained in [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF] and is equal to

                               k ( 15
)
where the values are expressed in SI units (N, m, rad). Let us simulate the calibration process assuming that the matrix (15) should be estimated by means of the identification algorithm described above, where the input data are generated by means of virtual experiments. In the frame of these experiments, the link is assumed to be fixed on one side and the external loading j W is applied on the another side. For each loading, the corresponding deflection vector is computed in accordance with expression • j j j  tW k ε , where j ε is the measurement noise. In accordance with the physical properties of the examined link and to conserve the linearity of the forcedeflection relation, the loading magnitude has been limited by 10N for the forces and 10Nm for the torques. The measurement noise magnitude has been defined as 25 for the orientation components (these values correspond to the precision of the best industrial measurement systems that currently are available on the market). These virtual experiments has been carried out six times, in order to obtain sufficient number of equations for the identification of 36 desired parameters ij k .

For these virtual experiments, the properties of the observation matrix used in the identification expression are quite good: rank is equal to 36 and the condition number is 1.00. Nevertheless, the identification are rather "surprising": the obtained compliance matrix essentially differs from the original one and is

                      k 6 6
                                (16)
Detailed comparison analysis of the original matrix k and its estimate k allows us to make the following conclusions concerning the harmful impact of the measurement noise on the identification of the elastostatic parameters in real industrial environment:

(i) the obtained compliance matrix k may lose the properties of positive-definiteness, which completely contradicts to the common physical sense that is based on the energybased definition of (iii) for some small elements, the identification accuracy may be extremely low (for example the element 11 k , which is 3 ~10 times less than 22 k and 33 k has been identified completely wrongly);

(iv) in the obtained matrix k , the number of non-zero elements is redundant compared to the original matrix k ; moreover, it is difficult to distinguish small elements ˆij k from socalled zero elements, which correspond to exact zeros in k induced by the physical properties of the examined link (for instance, the element 21 k that should be equal to zero by definition is the same order of magnitude as 11 k , which should be small but strictly positive);

(v) for the remaining elements, whose magnitude is high enough, the identification errors are quite acceptable (from 0.01% to 1.67%), but they should be further reduced by increasing number of the experiments.

It should be noted that for essentially lower measurement noise (with p  and   that are 100 times smaller) the above mentioned problems do not exist, however such measurement precision is not achievable in industrial environment at present.

Hence, as follows from this motivation example, the whole set of 36 elastostatic parameters   ij k composing the 6 6  matrix k cannot be estimated using commercially available measurement systems. The main reason for this difficulty is that, for some elements, corresponding deflections under the admissible loading are comparable with the measurement noise. To detect these indistinct elements, a simple indicator can be applied showing parameter-to-noise ratio (which is similar to signal-to-noise ratio in communication):

                      (17)
where ij  is a corresponding element of the relevant covariance

ij ij k                       (18)
allowing us to detect the same set of small or zero parameters whose identifiability is questionable. On the other side, the impact of these parameters on the elastostatic deflections is so small that they can be reasonably excluded from the desired stiffness model. These results confirm importance of the above pointed problems, which below are considered in details.

Research problem: practical identifiability

Summarizing theoretical background and simulation results presented in previous sub-sections, it is possible to make the following conclusions:

(i) complete elastostatic model of robotic manipulator includes huge number of parameters (258 for conventional 6 d.o.f. serial robot), whose simultaneous identification in presence of measurement noise is rather difficult or even impossible;

(ii) before applying the least-square identification technique , the manipulator elastostatic model should be reduced and redundant parameters should be eliminated, in order to ensure invertibility of the information matrix; this step can be performed using techniques similar to those developed for the geometrical calibration;

(iii) among the remaining non-redundant parameters, there are a number of non-significant ones, whose absolute values are relatively small, the identification accuracy is quite low and the impact on the compliance of the of the entire manipulator is almost negligible; these parameters can be treated as "practically non-identifiable" and should be also eliminated from the model, but relevant techniques are not available yet;

(iv) while developing relevant techniques allowing detection of "practically identifiable" parameters, it is prudent to take into account some specific properties of the compliance matrices induced by the elasticity physics such as the compliance matrix symmetry, presence of strictly zero elements (matrix sparseness), positive-definiteness, etc.

Hence, to obtain reliable stiffness model that is suitable for calibration, contains only significant and practically identifiable parameters while describing manipulator elastostatic properties sufficiently good, it is necessary to develop dedicated model reduction techniques and relevant rules allowing us to minimize number of parameters to be estimated and to reconstruct the original VJM-based model from these data taking into account mathematical relations between the model parameters caused by their physical sense.

PRACTICAL IDENTIFIABILITY IN MANIPULATOR CALIBRATION

Basic assumptions and terminology

Let us assume that the vector of desired elastostatic parameters π should be identified from the set of the linear equations [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF] whose least square solution is defined by the expression [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF], where the observation matrices ( , ) jj j  A q W are computed for certain set of measurement configurations   j q and loadings j W . Depending on the matrix set   j  A , corresponding system of linear equations can be solved for k π either uniquely or may have infinite number of solutions. In general, if the information matrix is rank-deficient, a general solution of the system (10) can be presented in the following form

  ˆ• •          πA I B A A λ . ( 19
)
where the superscript "+" denotes the Moore-Penrose pseudo inverse,

1 m TT j j j      A A A , 1 m TT j j j       B A
t and λ is an arbitrary vector of the same size as π . Using the later expression, all desired parameters contained in the vector π can be divided into the following groups [START_REF] Pashkevich | Computer-aided generation of complete irreducible models for robotic manipulators[END_REF]: G1: Identifiable parameters that can be obtained from [START_REF] Daney | Choosing measurement poses for robot calibration with the local convergence method and Tabu search[END_REF] in unique way and are independent from the arbitrary vector λ ;

G2: Non-identifiable parameters that cannot be computed uniquely from [START_REF] Daney | Choosing measurement poses for robot calibration with the local convergence method and Tabu search[END_REF] and can take on any value without influence on the right-hand side of the equation ( 9), they correspond to the zero columns of the observation matrix  A ;

G3: Semi-identifiable parameters that are also cannot be computed uniquely but have influence on the right-hand side of the equation ( 9); they are united in subgroups where a single one can be treated as identifiable if the remaining ones are fixed.

To obtain typical examples of the parameters belonging to the groups G1, G2 and G3, it is possible to use the ideas similar to geometrical calibration. For instance, the elastostatic parameters of the actuated joints and adjacent links are redundant in their totality and belong to the group G3. Besides, if the loading direction cannot be altered, a number of parameters belong to the group G2 and cannot be identified from the corresponding experimental data.

In this paper, in contrast to previous works, this classification is taking into account practical issues related to the limited precision of the measurement system. The main idea is to compare the absolute value of the estimated parameter with the range of possible fluctuations of the estimate caused by the measurement noise. For computational reasons, it is convenient to introduce a numerical indicator similar to the signal-to-noise ratio in communication, which is defined as follows

ˆ, 1, 2,... / i i i i     ( 20
)
where i  is the standard deviation of the parameter estimate ˆi  extracted from the diagonal of the covariance matrix (14). It is clear that i  can be treated as the inverse of the relative accuracy, which allows us to avoid the problem of division by zero. In the following sections this indicator will be referred to as parameter-to-noise ratio.

Using the above defined indicator, the set of parameters belonging to the group G1 (theoretically identifiable) can be further divided into three subgroups: G1+: Practically identifiable parameters, for which the accuracy indicator is high:

0 i  
 ; this subgroup describes principal elastostatic properties of the manipulator and should be certainly included in the reduced model used in the identification routines;

Copyright © 2013 by ASME

G1-:

Practically non-identifiable parameters, for which the accuracy indicator is low:

0 i  
 ; this subgroup contains non-essential parameters that can be assigned to zero in the VJM-model without essential impact on its precision (in practice, the majority of these parameters are nominally equal to zero due to the physical nature of the compliance matrices); G1~: Practically semi-identifiable parameters, for which the accuracy indicator is intermediate:

0 0 i     
; the parameters belonging to this subgroup are practically nonidentifiable for the current experimental setup but, hypothetically, can be converted into practically identifiable ones by increasing the experiment number, improving the measurement precision of by modification of the measurement configurations.

An open question however is related to justified assigning of the upper and lower bounds 0   and 0   . From practical point of view that is adopted below, it is reasonable to use 0 5

   and 0 2    , which is in a good agreement with the quantiles of the normal distribution. However, the user may modify these values in accordance with the specificity of the problem of interest.

The above presented definitions allow us to revise the concept of "suitable-for-calibration" model that in previous works included all parameters of the group G1 (this model is also referred to as the "complete and irreducible" one). In this work, this model is limited to include only parameters of the subgroup G1+ (practically identifiable) that can be estimated with reasonable accuracy and provide good approximation of the original complete model. The following subsections address different aspects of model reduction allowing us to obtain the desired model suitable for the elastostatic calibration.

It should be noted that, in spite of the fact that the main focus of the paper is on the elastostatic modelling, similar ideas can be also successfully applied in manipulator geometric calibration.

Model reduction: physical approach

Straightforward approach to the manipulator stiffness modeling leads to the exhaustive but redundant number of parameters to be identified. For instance, each links is described by a 6 6

 matrix that includes 36 parameters that are treated as independent ones. However, as follows from physics, number of the pure physical and independent parameters is essentially lower (for a trivial prismatic beam, for example, there are only five physical parameters: three describing the geometry and two describing the material properties). Hence, there are strong relations between these 36 parameters but this fact is usually ignored in elastostatic calibration. Besides, due to fundamental properties of conservative system, the desired compliance matrices should be strictly symmetrical and positive-definite. In addition, for typical manipulator links, the compliance matrices are sparsed due to the shape symmetry with respect to some axis, but this property is also not taken into account in identification of the elastostatic parameters.

To take advantages of the compliance matrix properties and to increase the identification accuracy, two simple methods can be applied that allows us to reduce the number of parameters to be computed in the identification procedure [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. They can be treated as the physics-based model reduction techniques and formalized in the following way. M1:Symmetrisation. For all compliance matrices k to be identified, replace the pairs of symmetrical parameters

  , ij ji
kk by a single one ,

ij k i j  .
For each link, this reduction procedure is equivalent to redefinition of the model parameters vector in the following way

•   π M π (21)
where the binary matrix M of size 36 21  describes the mapping from the original to reduced parameter space. It can be proved that corresponding basic expression for the identification (9) can be rewritten as

)• (,      A q W t π (22) where (.) (.

)•    A

A M denotes the reduced observation matrix. The later can be also computed as

θ 1 θ θ 2 θ θ 21 θ ( , ) , , ... [ ] T T T    A q W J ω J W J ω J W J ω J W (23)
where 12 , ,... ωω denote the binary matrices of size 6 6  for which non-zero elements (i.e. equal to 1) are located in the following way: for the parameter l  corresponding to the matrix elements , ij k i j  , the non-zero elements are 1 ij ji   . It is clear that this idea allows us to reduce the number of links compliance parameters from 36 to 21 (and from 258 to 153 for the entire 6 d.o.f. manipulator).

M2:

Sparcing. For all compliance matrices k to be identified, eliminate from the set of unknowns the parameters ij k corresponding to zeros in the stiffness matrix template 0 k derived analytically for the manipulator link with similar shape.

To derive the desired template matrix it is convenient to use any realistic link-shape approximation. For example using the trivial beam [START_REF] Timoshenko | Theory of elasticity[END_REF], the desired template can be presented as 0 *00000 0*000* 00*0*0 000*00 00*0*0 0*000*

         k ( 24
)
where the symbol "*" denotes non-zero elements. It allows further reducing the number of the unknown parameters from 21 to 8, taking into account only essential ones from physical point of view. It can be also proved that the template (24) is valid for any link whose geometrical shape is symmetrical with respect to three orthogonal axes. But it is necessary to be careful if this property is not kept strictly.

Summarizing these methods, it should be mentioned that the above presented methods essentially reduce the number of parameters to be identified (by the factor 4.5) but they do not violate such basic properties as the mode completeness, i.e. the ability to describe any deflection caused by the external loading. However, the reduced model may still have some redundancy in the frame of entire manipulator, where the virtual springs of adjacent joints/actuators cause similar impact on the end-effector deflections under the loading. To eliminate these parameters (belonging to the group G3), relevant algebraic technique developed in our previous work [START_REF] Pashkevich | Computer-aided generation of complete irreducible models for robotic manipulators[END_REF] can be applied.

To illustrate efficiency of the proposed methods, the identification problem considered in section 3.2 have been solved for reduced set of the compliance parameters. It yielded the following result

ˆ•                                 k (25)
which is essentially better compared to [START_REF] Pashkevich | Computer-aided generation of complete irreducible models for robotic manipulators[END_REF]. In particular, the identification errors for the most of the desired parameters are less than 0.4%, i.e. 4 times lower. The only exception is the small element 11 k that is still negative and contradicts to the physical sense. This motivates further efforts to obtain reliable model suitable for elastostatic calibration in industrial environment.

Model reduction: statistical approach

As follows from relevant study and above presented example, rigorous reduction methods based on the physical and mathematical properties of the compliance matrix are rather limited if the measurement noise is non-negligible. This gives us reasons to develop some heuristic rules that take into account the measurement noise impact on the identification accuracy. It is clear that extremely low accuracy is not acceptable, but often corresponding parameters are so small that their influence on the end-effector deflections is almost negligible. This supports an idea for heuristic reduction of small model parameters but leaving an open problem of their further reconstruction in the VJM-model using some empirical or semi-empirical relations induced by mathematical relations between the stiffness matrix elements.

To take into account the relative accuracy of the parameter estimates, it is convenient to use a simple indicator showing parameter-to-noise ratio [START_REF] Timoshenko | Theory of elasticity[END_REF] introduced in sib-section 4.1. It is evident that it should be applied only to those parameters that belong to the group G1 (theoretically identifiable). Using this index, a heuristic model reduction technique allowing us to distinguish the practically identifiable parameters from the hardly-identifiable ones can be formalized as follows: M3:Neglecting.

Step 1:

Using complete but non-redundant model derived after application of physical and algebraic model reduction techniques, compute estimates of the desired parameters ˆ π and their covariance matrix   ĉov  π by means of equations ( 12) and ( 14);

Step 2: Using the parameters estimates ˆ π and the diagonal elements of the covariance matrix   ĉov  π , compute the parameter-to-noise ratios i  in accordance with expression (20);

Step 3. For all compliance matrices k to be identified, eliminate from the set of unknowns the parameters ij k for which parameter-to-noise ratios ij  is lower the user defined threshold:

0 ij    .
This method allows us to eliminate from the model the parameters whose identification accuracy is comparable with the noise impact and, strictly speaking, these values can not be considered as reliable estimates of ij k .

As follows from our experience, it a very powerful method with two useful features: (i) elimination of small (but theoretically non-zero) parameters, and (ii) detection of elements corresponding to zeros in the matrix template (see method M2), if the latter has been defined rather carefully. These conclusions are clearly confirmed by the numerical example presented in section 3.2 (see [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods, and Algorithms[END_REF] and ( 18)). In this example, it is worth to pay attention to the element 11 k that is really small for the majority of manipulator links (because the link length is always essentially higher compared to the crosssection dimensions). Elimination of this parameter is really negligible for the manipulator Cartesian stiffness matrix (7) that integrates impact of all compliance elements. Nevertheless, after identification, the parameter 11 k can be reconstructed approximately using non-zero elements of the compliance matrix and some relations between ˆij k induced by physics. The last problem is currently under study but is not in the scope of this paper.

In conclusion of this section, it should be noted that the proposed methods M1-M3 allow us essentially reducing the number of model parameters while retaining the model accuracy. For example, for 6 d.o.f. manipulator, the number of Copyright © 2013 by ASME parameters is reduced from 258 to 42 allowing us to obtain an adequate stiffness model in real industrial environment.

APPLICATION EXAMPLE

The developed model reduction techniques have been applied to the elastostatic identification of the industrial robot KUKA KR-270 (Figure 3), which is used for processing of aircraft components made from contemporary high performance materials. The application area requires high positional precision while the technological process generates significant end-effector deflections, which can be compensated on the control level via the elastostatic modeling. It should be noted that for this robot, the elastostatic model must also take into account the impact of the gravity compensator, which creates the closed loop between the first and second links and is described in the stiffness model by two additional parameters (the spring stiffness and preloading). The geometrical parameters of the gravity compensator have been identified from dedicated experiments. For the elastostatic calibration, two different experiment setups have been used. The first of them is based on the application of the gravity force 2.5 kN to the special measurement tool, which contains three reference points (Figure 4). Relevant deflections have been measured by the laser-tracer Leica (absolute accuracy 10μm). For the second setup, the deflections have been generated by applying horizontal loadings 3.5 kN. To find the measurement configurations that ensure the best identification accuracy, the design of experiments technique has been applied, which is based on the dedicated industry-oriented performance measure proposed in our previous work [START_REF] Pashkevich | Computer-aided generation of complete irreducible models for robotic manipulators[END_REF]. These optimal configurations have been obtained taking into account physical constraints that are related to the joint limits and the work-cell obstacles. For the considered manipulator, the complete elastostatic model includes 260 parameters (with 2 additional parameters for the gravity compensator). To obtain the suitable for calibration model, the original set of parameters has been sequentially reduced by applied physical, algebraic and statistical methods presented in Section 4. This yielded the simplified model with 18 parameters that can be identified in practice with the desired accuracy (the most essential of them corresponding to the joint actuators are presented in Table 1). Relevant model has been further used for the compliance error compensation for the robotic based milling of a aircraft parts, where essential improvement of the precision has been achieved. For this application, the reduced elastostatic model allowed us to compensate more than 95% of elastostatic deflections and to ensure positional accuracy about 0.1 mm (that is comparable with the robot repeatability 0.06 mm).

CONCLUSIONS

The paper deals with the problem of the manipulator stiffness modeling, which is extremely important for the robotic-based machining of contemporary aeronautic materials where high position accuracy is required while performing prescribed manufacturing task. The main attention is paid to the elastostatic parameters identification and model reduction, where the notion of practical identifiability is introduced that relies on the essential differences in the model parameter magnitudes and the measurement noise impact.

In contrast to previous works, the manipulator stiffness properties are described by the sophisticated model, which takes into account the flexibilities of all mechanical elements such as links, actuated joints, mechanical transmissions, etc.. In the frame of this model, the virtual joint method (VJM) is used, which operates with 6×6 stiffness matrices for each compliant link and scalar coefficients for the joints/transmissions. This

Copyright © 2013 by ASME yields extremely high number of the model elastostatic parameters to be identified, that for a conventional 6 d.o.f. manipulator reaches 258. Even by eliminating dependent ones does not allow us to reduce this number substantially; for this manipulator the stiffness model includes 153 independent parameters that, theoretically, may be identified. However, the parameter magnitudes differ significantly (~1000 times), so straightforward application of conventional identification technique does not give reliable results (for some parameters the estimation errors are greater than 100% that also may violate fundamental physical properties of the stiffness matrices, such as positive-definiteness and symmetry). On the other hand, some of the desired parameters are so small that their influence on the manipulator accuracy is negligible. This leads the problem of further reduction of the stiffness model that aims at eliminating some small parameters. To distinguish these small parameters from essential ones, the notion of practical identifiability is introduced.

To solve the problem, physical and statistical model reduction methods are developed. They take into account mathematical relations between the elements of the compliance matrices and parameter magnitude with respect to the measurement noise impact. The advantages of the developed approach are illustrated by an application example that deals with the stiffness modeling of industrial robot used in aerospace industry. In future, the problem of the complete model reconstruction from the obtained set of practically identifiable parameters will be in the focus of our attention.

Figure 2

 2 Figure 2 Manipulator link considered in the motivation example

Figure 3 Figure 4

 34 Figure 3 Experimental setup for the identification of the elastostatic parameters x

Table 1

 1 Principal elastostatic parameters of KR-270

	Parameter	Value	Accuracy
	k c , [rad•μm/N]	0.144 ±0.031	21.5%
	s 0 , [mm]	458 ±27	5.9%
	k 2 , [rad•μm/N]	0.302 ± .004	1.3%
	k 3 , [rad•μm/N]	0.406 ±0.008	2.0%
	k 4 , [rad•μm/N]	3.002 ±0.115	3.8%
	k 5 , [rad•μm/N]	3.303 ±0.162	4.9%
	k 6 , [rad•μm/N]	2.365 ±0.095	4.0%

ACKNOWLEDGMENTS

The work presented in this paper was partially funded by the Region "Pays de la Loire", France and by the projects ANR COROUSSO and FEDER ROBOTEX, France.


Alexandr Klimchik
email: alexandr.klimchik@mines-nantes.fr

Stphane Caro
email: stephane.caro@irccyn.ec-nantes.fr

Yier Wu
email: yier.wu@mines-nantes.fr

Damien Chablat
email: damien.chablat@irccyn.ec-nantes.fr

Benoit Furet
email: benoit.furet@univ-nantes.fr

Anatol Pashkevich
email: anatol.pashkevich@mines-nantes.fr

Stiffness modeling of robotic manipulator with gravity compensator

Keywords: Stiffness modeling, gravity compensator, industrial robot

The paper focuses on the stiffness modeling of robotic manipulators with gravity compensators. The main attention is paid to the development of the stiffness model of a spring-based compensator located between sequential links of a serial structure. The derived model allows us to describe the compensator as an equivalent non-linear virtual spring integrated in the corresponding actuated joint. The obtained results have been efficiently applied to the stiffness modeling of a heavy industrial robot of the Kuka family.

Introduction

in its turn, influences on the position accuracy and stiffness properties. The latter motivates enhancement of the manipulator stiffness modeling techniques that are in the focus of this paper. The problem of stiffness modeling for serial manipulators has been studied in the robotic community from different aspects [START_REF] Company | A Method for Modeling Analytical Stiffness of a Lower Mobility Parallel Manipulator[END_REF][START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF][START_REF] Tyapin | Kinematic and elastostatic design optimization of the 3-DOF Gantry-Tau parallel kinamatic manipulator[END_REF]. In particular, special attention has been paid to heavy industrial and medical robots with essential deflections of the end-effector [START_REF] Yi | Geometric analysis antagonistic stiffness redundantly actuated parallel mechanism[END_REF][START_REF] Meggiolaro | Geometric and elastic error calibration of a high accuracy patient positioning system[END_REF][START_REF] Kvecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF]. Among a number of existing stiffness modeling ap-proaches, the Virtual Joint Modeling (VJM) method looks the most attractive for the robotic applications. Its main idea is to take into account the elastostatic prop-erties of flexible components by presenting them as equivalent virtual springs lo-calized in the actuated or passive joints [START_REF] Salisbury | Active stiffness control of a manipulator in Cartesian coordinates[END_REF]. Because of its simplicity and effi-ciency, this approach has been applied to numerous case-studies and progressively enhanced to take into account specific features of robotic manipulators, such as in-ternal and external loadings, closed-loops, etc. [START_REF] Klimchik | Compensation of tool deflection in robotic-based Milling[END_REF][START_REF] Quennouelle | Instantaneous Kinemato-Static Model of Planar Compliant Parallel Mechanisms[END_REF]. Besides, some authors have extended the VJM approach by using advanced 6 d.o.f. virtual springs describing stiffness properties of the manipulator elements [START_REF] Pashkevich | Accuracy Improvement for Stiffness Modeling of Parallel Manipulators[END_REF]. However, the problem of the stiffness modeling of the manipulators with gravity compensators has not been studied yet sufficiently; there are only limited number of works that addressed this issue. Besides, in previous works, the main attention has been paid to the compen-sator design [START_REF] Takesue | Design and Prototype of Variable Gravity Compensation Mechanism (VGCM)[END_REF] and modification of the inverse kinematics algorithms integrated in the robot controller [START_REF] Luca | A PD-type Regulator with Exact Gravity Cancellation for Robots with Flexible Joints[END_REF]. In contrast to previous works, this paper deals with the stiffness modeling of serial robotic manipulators equipped with a springbased compensator located be-tween sequential links. It proposes an extension of the VJM-based approach al-lowing to integrate the elasto-static properties of the gravity compensator (that creates a closed loop) into the conventional stiffness model of the manipulator.

Stiffness model of gravity compensator

Let us consider the general type of gravity compensator that incorporates a passive spring attached between two sequential links of the manipulator (Fig. 1). In a such architecture, there is a closed loop that generates an additional torque in the manipulator joint. So, the specificity of such design allows us to limit modifications of the robot stiffness model by adjusting the virtual joint stiffness parameters.

The geometrical model of considered compensator includes three main node points P 0 , P 1 , P 2 . Let us denote corresponding distances as L = |P 1 , P 2 |, a = |P 0 , P 2 |, s = |P 0 , P 1 |. In addition, let us introduce the angles α, ϕ and the distances a x and a y , whose geometrical meaning is clear from Fig. 1. Using these notations, the compensator spring deflection s can be computed as

s = a 2 + L 2 + 2 • a • L • cos(α + q) (1)
and evidently depends on the joint angle q. The considered mechanical design allows us to balance the manipulator weight for any given configuration by adjusting the compensator spring preloading. It can be taken into account by introducing the zero-value of the compensator length s 0 corresponding to the unloaded spring. Under this assumption, the compensator force applied to the node P1 can be expressed as

F s = K c • (s -s 0 )
, where K c is the compensator spring stiffness. Further, the angle ϕ between the compensator links P0P1 and P1P2 (see Fig. 1) can be found from the expression sin ϕ = a/s • sin(α + q), which allows us to compute the compensator torque M c applied to the second joint as

M c = K c • (1 -s 0 /s) • a • L • sin(α + q) (2)
After differentiation of the latter expression with respect to q, the equivalent stiffness of the second joint (comprising both the manipulator and compensator stiffnesses) can be expressed as:

K θ = K 0 θ + K c • a L • η q (3)
where the coefficient

η q = cos(α + q) - s 0 s • a L s 2 sin 2 (α + q) + cos(α + q) (4)
highly depends on the value of joint variable q and the initial preloading in the compensator spring described by s 0 . Hence, using expression (3), it is possible to extend the classical stiffness model of the serial manipulator by modifying the virtual spring parameters in accordance with the compensator properties. While in the paper this approach has been used for the particular compensator type, similar idea can be evidently ap-plied to other compensators. It is also worth mentioning that the geometrical and elastostatic parameters of gravity compensators (α, a, L and K c , s 0 for the presented case) usually are not included in datasheets. For this reason, these parameters should be identified from the calibration experiments (see Section IV).

Extension of the VJM-based approach

Stiffness of a serial robot highly depends on its configuration and is defined by the Cartesian stiffness matrix. Using the VJM-based approach adopted in this paper, the manipulator can be presented as the sequence of rigid links separated by the actuators and virtual flexible joints incorporating all elastostatic properties of flexible elements [START_REF] Kvecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF]. In accordance with these assumptions, the elastostatic model of the serial robot can be presented as shown in Fig. For this serial robot, the end-effector location t can be defined by the vector function g (q, θ ), where q, θ denote the vectors of the actuator and virtual joint coordinates respectively. It can be proved that the static equilibrium equations can be written as J T θ • F = K θ • θ , where J θ = ∂ g (q, θ )/∂ θ is the Jacobian matrix, the matrix K θ = diag K θ 1 , ..., K θ n aggregates stiffnesses of all virtual springs and F is the external loading applied to the robot end-effector. In order to find the desired stiffness matrix K C , the force-deflection relation should be linearized in the neighborhood of the current configuration q and presented in the form F = K C • ∆ t, where ∆ t is the end-effector deflection caused by the external loading F. After relevant transformations, one can get the following expression for the desired stiffness matrix [START_REF] Yi | Geometric analysis antagonistic stiffness redundantly actuated parallel mechanism[END_REF]

K C = J θ • K -1 θ • J T θ -1 (5)
which depends on the manipulator geometry and its elastostatic properties. This classical expression has been originally derived for the case of manipulators without gravity compensators, where the matrix K θ is constant. To integrate into this model the gravity compensator, it is proposed to replace the classical joint stiffness matrix K 0 θ (that takes into account elastostatic properties of serial manipulator without gravity compensator) by the sum of two matrices

K θ = K 0 θ + K GC θ (q) (6)
where the second term K GC θ (q) is configuration-dependent and takes into account the elastostatic properties of the gravity compensator (by transforming them in the virtual joint space). In practice, the matrix K GC θ (q) may include either one or several non-zero elements (depending on the number of compensators). These elements ale located on the matrix diagonal and are computed using [START_REF] Tyapin | Kinematic and elastostatic design optimization of the 3-DOF Gantry-Tau parallel kinamatic manipulator[END_REF]. Hence, by modification of the manipulator joint stiffness matrix in accordance with expression [START_REF] Kvecses | The stiffness matrix in elastically articulated rigid-body systems[END_REF], it is possible to integrate the compensator parameters into the classical VJM-based model of the serial manipulator. It should be mentioned that this idea can be also applied for other types of compensators.

Application example

Let us illustrate the efficiency of the developed stiffness modeling technique by applying it to the compliance error compensation in the robotic-based machining performed by the industrial robot KUKA KR-270. This robot is equipped with the spring-based gravity compensator located between the first and second links (which leads to the influence on the second actuated joint). In accordance with the considered specifications, the technological process should be performed in the square area of the size 2000 mm2000 mm located at the height 500 mm over the floor level (see Fig, 3 for more details). For comparison purpose, it is assumed that machining force is constant throughout the working area and it is equal to F = (0, 360N, 560N, 0, 0, 0) T , which corresponds to a typical milling process. To obtain the desired stiffness model, special calibration experiments have been conducted and dedicated identification procedures have been applied. These yield the desired geometrical and elastostatic parameters of the robot with gravity com-pensator presented in Table 1. From these data, the classical joint stiffness matrix K 0 θ has been constructed straightforwardly (from parameters k 1 , ...k 6 ). To integrate the gravity compensator properties, the matrix K GC θ (q) has been computed using expression (3) and values k c , L, a x and a y for any robot configuration q. Relevant computations have been done throughout the required working area, where the endeffector deflections has been evaluated for given machining force F. The computational results are summarized in Fig. 4, where the end-effector compliance errors are presented. These results show that the compliance errors are not negligible here and vary from 0.34 to 3.5 mm. So, to achieve the desired precision, it is reasonable to apply the error compensation technique, which is based on the reliable stiffness model. For comparison purposes, two alternative stiffness models have been examined. The first one is based on the classical assumptions and takes into account the stiffness properties of the actuated joints only (without gravity compensator). The second model integrates the gravity compensator in accordance with the approach proposed in this paper. Both the models have been used for the compliance error compensation for the described above manufacturing task. Relevant results are presented in Fig. 5, which shows the difference in the compliance error compensation while applying the classical and the proposed approaches. As follows from the figure, ignoring the gravity compensator influence may lead to the position error over/under compensation of the order of 0.07 mm, which is not admissible for the manufacturing processes employed in the aerospace applications studied here. Hence, to ensure the high precision for the robotic-based machining, the compliance error compensation technique must rely on the stiffness model, which takes into account the impact of the gravity compensator, in accordance with the approach developed in this paper.

Conclusions

The paper presents a new approach for the stiffness modeling of robots with the spring-based gravity compensators, which are located between the manipulator se-quential links. Using this approach, the compensator has been replaced by an equivalent non-linear virtual spring integrated in the corresponding actuated joint. This methodology allowed us to extend the VJM-based modeling technique for the case of manipulators with closed-loops induced by the gravity compensators via using configuration-dependent joint stiffness matrix. Efficiency of the developed approach and its industrial value have been confirmed by an application ex-ample, which deals with robotic-based milling of large-dimensional parts for aerospace industry. Feature work will deal with integration of this modeling approach into the robotic software CAD system.

Fig. 1

 1 Fig. 1 Mechanics of gravity compensator.

2

 2

Fig. 2

 2 Fig. 2 Serial industrial robot (a) and its VJM-based stiffness model (b).

Fig. 3

 3 Fig. 3 Industrial robot Kuka KR270 with required machining area.

Fig. 4

 4 Fig.4The end-effector deflections caused by machining forces throughout working area.

Fig. 5

 5 Fig. 5 Difference in stiffness error compensation between two strategies.

Table 1

 1 Elastostatic and geometrical parameters of robot with gravity compensator.

	Parameters	Units	Values	CIs
	k 1 k 2 k 3 k 4 k 5 k 6 k c	[rad • m/N] [rad • m/N] [rad • m/N] [rad • m/N] [rad • m/N] [rad • m/N] [rad • m/N]	3.774 • 10 -6 0.302 • 10 -6 0.406 • 10 -6 3.002 • 10 -6 3.303 • 10 -6 2.365 • 10 -6 0.144 • 10 -6	-0.004 • 10 -6 0.008 • 10 -6 0.115 • 10 -6 0.162 • 10 -6 0.095 • 10 -6 0.031 • 10 -6
	s 0	[mm]	458.00	27.0
	L	[mm]	184.72	0.06
	a x	[mm]	685.93	0.70
	a y	[mm]	120.30	0.69

Acknowledgements The work presented in this paper was partially funded by the ANR, France (Project ANR-2010-SEGI-003-02-COROUSSO) and Project ANR ROBOTEX. The authors also thank Fabien Truchet, Guillaume Gallot, Joachim Marais and Sbastien Garnier for their great help with the experiments.


Keywords: integer quadratic knapsack problem, separable objective function, direct expansion, binary expansion, piecewise interpolation

linear interpolation can not be applied to transform (QMKP) into an equivalent 0-1 piecewise linear problem.

Rewriting integer variables into zero-one variables: some guidelines for the integer quadratic multi-knapsack problem

Abstract: This paper is concerned with the integer quadratic multidimensional knapsack problem (QMKP) where the objective function is separable. Our objective is to determine which expansion technique of the integer variables is the most appropriate to solve (QMKP) to optimality using the upper bound method proposed by [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. To the best of our knowledge the upper bound method previously mentioned is the most effective method in the literature concerning (QMKP). This bound is computed by transforming the initial quadratic problem into a 0-1 equivalent piecewise linear formulation and then by establishing the surrogate problem associated. The linearization method consists in using a direct expansion initially suggested by [START_REF] Glover | Improved linear programming formulations of nonlinear integer problems[END_REF] of the integer variables and in applying a piecewise interpolation to the separable objective function. As the direct expansion results in an increase of the size of the problem, other expansions techniques may be utilized to reduce the number of 0-1 variables so as to make easier the solution to the linearized problem. We will compare theoretically the use in the upper bound process of the direct expansion (I) employed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] with two other basic expansions, namely: (II) a direct expansion with additional constraints and (III) a binary expansion. We show that expansion (II) provides a bound which value is equal to the one computed by [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Conversely, we provide the proof of the non applicability of expansion (III) in the upper bound method. More specifically, we will show that if (III) is used to rewrite the integer variables into 0-1 variables then a

Introduction

This paper deals with the integer quadratic multi-knapsack problem (QMKP) where the objective function is separable. Problems of this structure arise in numerous industrial and economic situations, for instance in production planning [START_REF] Stecke | Nonlinear integer production planning problems[END_REF], reliability allocation [START_REF] Nakagawa | Optimal reliability allocation by branch and bound techniques[END_REF] and finance [START_REF] Faaland | An integer programming algorithm for portfolio selection[END_REF]. These include the main application of (QMKP) which is in the portfolio management area where the investments are independent, see [START_REF] Djerdjour | A surrogate-based algorithm for the general quadratic multidimensional knapsack[END_REF] and [START_REF] Faaland | An integer programming algorithm for portfolio selection[END_REF]. Nevertheless, solving (QMKP) efficiently will constitute a starting point to solve the more general and realistic portfolio management problem where the investments are dependent, i.e. the objective function is non separable.

The integer quadratic multi-knapsack problem (QMKP) where the objective function is separable consists in maximizing a concave separable quadratic integer function subject to m linear capacity constraints. It may be stated mathematically as follows:

( ) ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ = = ≤ ≤ = ≤ = - = ∑ ∑∑ = = = n j x n j u x m i b x a x f x d x c x f QMKP j j j n j i j ij n j n j j j j j j j , , 1 integer, , , 1 , 0 , , 1 , t. s.

) ( max ) (

where the coefficients c j , d j , a ij , b i are nonnegative. The bounds u j of variables x j are pure integers, with ( )

j j j d c u 2 / ≤
. Indeed, the separable objective function is concave which implies that for all function f j , ( )

j j j d c x 2 / * ≤
, where x * j is the optimal solution of the program ( )

j j x x f j 0 max ≥ .
The problem (QMKP) which is a NP-hard problem [START_REF] Chaillou | Best network flow bounds for the quadratic knapsack problem[END_REF] is a generalization of both the integer quadratic knapsack problem [START_REF] Bretthauer | The nonlinear knapsack problem -algorithms and applications[END_REF] and the 0-1 quadratic knapsack problem where the objective function is subject to only one constraint [START_REF] Billionnet | A new upper bound for the 0-1 quadratic knapsack problem[END_REF].

Since, (QMKP) is NP-hard, one should not expect to find a polynomial time algorithm for solving it exactly. Hence, we are usually interested in developing branch-andbound algorithms. A key step in designing an effective exact solution method for such a maximization problem is to establish a tight upper bound on the optimal value. Basically, the available upper bound procedures for (QMKP) may be classified as attempting either to solve efficiently the LP-relaxation of (QMKP) (see [START_REF] Bretthauer | The nonlinear knapsack problem -algorithms and applications[END_REF] and [START_REF] Körner | Integer quadratic programming[END_REF]) or to find a good upper bound, of better quality than the LP-relaxation of (QMKP), transforming (QMKP) into a 0-1 linear problem easier to solve (see [START_REF] Djerdjour | A surrogate-based algorithm for the general quadratic multidimensional knapsack[END_REF] and [START_REF] Mathur | A branch and search algorithm for a class of nonlinear knapsack problems[END_REF]). To the best of our knowledge, the upper bound method we have proposed in a previous work [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] is better than the existing methods (Djerdjour, Mathur and Salkin algorithm [START_REF] Djerdjour | A surrogate-based algorithm for the general quadratic multidimensional knapsack[END_REF], a 0-1 linearization method, a classical LP-relaxation of (QMKP)) from both a qualitative and a computational standpoint. We have first used a direct expansion of the integer variables, originally suggested by Glover [START_REF] Glover | Improved linear programming formulations of nonlinear integer problems[END_REF], and apply a piecewise interpolation to the objective function: an equivalent 0-1 linear problem is thus obtained. The second step of the algorithm consists in establishing and solving the surrogate relaxation problem associated to the equivalent linearized formulation.

Nevertheless the transformed linear formulation encounters numerous 0-1 variables because of the direct expansion used (denoted by expansion (I) in the following).

Consequently, other expansions techniques may be utilized to reduce the number of 0-1 variables so as to make easier the solution to the linearized problem. Let us consider the three basic expansions for rewriting the integer variables of (QMKP) into 0-1 variables:

• Expansion (I): direct expansion

∑ = = j u k j jk x y 1 , { } j j u x K 1 ∈ , { } 1 , 0 ∈ jk y n j K 1 = ∀ , j u k K 1 = • Expansion (II): direct expansion with additional constraint ∑ = = j u k j jk x y k 1 ' , ∑ = ≤ j u k jk y 1 ' 1 { } j j u x K 1 ∈ , { } 1 , 0 ' ∈ jk y n j K 1 = ∀ , j u k K 1 = • Expansion (III): binary expansion. ⎡ ⎤ ∑ = = ) log( 1 2 j u k j jk k x z , ⎡ ⎤ ∑ = ≤ ) log( 1 2 j u k j jk k u z { } j j u x K 1 ∈ , { } 1 , 0 ∈ jk zn j K 1 = ∀ , ⎡ ⎤ ) log( 1 j u k K =
The purpose of this note is to evaluate the impact of the use of the above expansion techniques, on the upper bound computation developed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. More specifically, we will determinate which expansion is the most appropriate to be used in the upper bound method for (QMKP) [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. We will compare theoretically the use of the direct expansion (I) with the direct expansion with additional constraints (II) and with the binary expansion (III). We will show that the use of (II) provides a bound which value is equal to the one computed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Conversely, we provide the proof of the non applicability of both (III) and a linear interpolation to transform (QMKP) into an equivalent 0-1 piecewise linear problem.

The paper is organized as follows. The next section summarizes the upper bound method developed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] detailing the direct expansion (I) of integer variables and the piecewise interpolation. In Section 3, the direct expansion with additional constraints (II) is applied to (QMKP) so as to compute the upper bound suggested by Quadri et al. [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Section 4 is dedicated to the binary expansion (III). We finally conclude in Section 5.

In the remainder of this paper, we adopt the following notations: letting (P) be an integer or a 0-1 program, we will denote by ( P ) the continuous relaxation problem of (P). We let Z[P] be the optimal value of the problem (P) and Z[ P ] the optimal value of ( P ). Finally ⎡⎤ x (resp. ⎣⎦ x ) will denote the smallest (resp. highest) integer greater (resp. lower) than or equal to x.

Section 2. Direct expansion of the integer variables

In this section we summarize the upper bound method for (QMKP) proposed by Quadri et al. [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. First, an equivalent formulation is obtained by using a direct expansion (I) of the integer variables x j as originally proposed by Glover [START_REF] Glover | Improved linear programming formulations of nonlinear integer problems[END_REF] and by applying a piecewise interpolation to the initial objective function as discussed in [START_REF] Djerdjour | A surrogate-based algorithm for the general quadratic multidimensional knapsack[END_REF].

The direct expansion of the integer variables x j consists in replacing each variables x j by a sum of u j 0-1 variables y jk such that ∑ = = In the second step of the algorithm, a surrogate relaxation is applied to the LPrelaxation of (MKP). The resultant formulation (KP, w) is the surrogate relaxation problem of (MKP) and can be written as: As proved by Glover [START_REF] Glover | Improved linear programming formulations of nonlinear integer problems[END_REF], (KP, w) is a relaxation of (MKP). For any value of

( ) ( ) {} ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = = ∈ = ≤ = ∑∑ ∑∑ == == 1 1 , 1 , 0 , , 1 , t. s.
( ) [ ] {} ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = = ∈ ≤ ∑∑ ∑ ∑ ∑∑ == = = == 1 1 , 1 , 0 t. s. max ) , (
0 ≥ w [ ] w KP Z ,
the optimal value of ( ) proposed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF].

Section 3. Direct expansion of integer variables with additional constraints

In this section we apply to the integer variables of (QMKP) the direct expansion with additional constraints (II). That is each variable x j is replaced by the following expression ∑ = ( ) Following the upper bound method developed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] we then establish the surrogate problem (KP 2 ,w) associated to (MKP 2 ). The problem (KP 2 ,w) can be written as: Utilizing the optimal solution * ω of the dual of ) ( 2MKP as the surrogate multiplier, the optimal value of the LP-relaxation of (KP 2 , * ω ) provides an upper bound of (QMKP). The following proposition and its corollary show that this upper bound, computed in this section through the use of expansion (II), is equal to the one computed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. . We denote:

( ) {} ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ = = = ∈ ≤ = ≤ = ∑ ∑∑ ∑∑ = == == 1 1 ,..., 1 , 1 , 0 1 , , 1 , t. s.
-

∑ = = j u k jk j y k x 1 ' . Obviously [ ] ; , 0 j j u x ∈ - ⎣ ⎦ j j x p =
, where ⎣ ⎦

x stands for the highest integer lower than or equal to the real x. j p is the integer part of j x . Note that j j u p < since . This implies that y is feasible for ( ) w KP, since we have: ))

∑ ∑∑ ∑ ∑∑ ∑∑ ∑ = == = == == = ≤ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ m i i i
{ } ∑ ∑ = = ≤ ∈ ∀ j j u k jk jk u k jk jk y s y f n j 1 1 ' , , , 1 K Let { } n j , , 1 K ∈
. We first note that:

( ) ( ) ( ) 1 , , , 1 , 1 1 , 1 , 1 1 1 + + = - + = = + - = - + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ - = + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ = ∑ ∑ ∑
K ε ( ) ( ) ⎪ ⎩ ⎪ ⎨ ⎧ ≥ = ≥ + + + 0 , , 1 t. s. minimize α β α β ε α j jk j j j u k f k p D K 1 0 1 t. s. ) ' ( maximize ' 1 ' 1 ' 1 '
The solution y ~ defined by (

{ }, , , 0 j u k K ∈ ∀ , jk f k ≥ + β α ~ Let us denote ( ) jk p j j p j j j f f p k f k p k j j - - + - + = Π +1 , , ) ( 1 ) (
. We thus have to prove that:

{ }, , , 0 j u k K ∈ ∀ , 0 ) ( ≥ Π k j Using the definition of jk f ( 2 k d k c f j j jk - =
), we find with some easy algebra:

) 1 ( ) 1 2 ( ) ( 2 + + + - = Π j j j j j j j p p d k d p k d k Thus, ) (k j Π
is a degree 2 polynomial in k. Its discriminant is equal to Let us now prove that ) ' ( ) ( y h y g = .

{ }, , , 0 j u k K ∈ ∀ 0 ) ( ≥ Π k j , which implies that ( )
-For 0 j j = , we have: We thus have proved in this Section that solving (QMKP) with a linearization technique and a surrogate relaxation may be equivalently done with expansion (I) or with expansion (II) of the integer variables. We proved that expansions (I) and (II) are equivalent in the sense that they provide the same upper bound. However, expansion (II) involves n added constraints without improving the quality of the upper bound in comparison with expansion (I).

( ) 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 , , ' , 1 , ' , , 1 ' , ,

Section 4. Binary expansion of integer variables

This section is dedicated to the use of a binary expansion (referred as expansion (III))

of this integer variables in the upper bound procedure developed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Such expansion consists in rewriting each integer variable x j (j=1,…,n) as: A possible way to get a further decrease of the number of 0-1 variables would be the use of a classical linearization technique as suggested by [START_REF] Foret | Applications de l'algèbre de Boole en recherché opérationnelle[END_REF] [START_REF] Foret | Applications de l'algèbre de Boole en recherché opérationnelle[END_REF]. This linearization scheme would first consist in applying (III) to the integer variables and then using the basic linearization technique for 0-1 variables. Nevertheless, preliminary computational experiments have shown that the bound provided by this technique is of worst quality than the one computed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Moreover, this process seems to be more time consuming than [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF].

(

  LAMSADE, Université Paris Dauphine, Place du Maréchal de Lattre de Tassigny, 75775 Paris Cedex 16, France) Eric Soutif (CEDRIC, Conservatoire National des Arts et Métiers, 292 rue Saint Martin, 75141 Paris Cedex 03, France)

1 .

 1 Since the objective function f is separable a linear interpolation can then be applied to each objective function term f j . Consequently, (QMKP) is equivalent to the 0-1 piecewise linear program (MKP):

  w KP, constitutes an upper bound for [ ] MKP Z . But the value of the bound [ ] w KP Z ,depends on the choice of the surrogate multiplier employed. It is proved in[START_REF] Nakagawa | Optimal reliability allocation by branch and bound techniques[END_REF] that if w* is chosen as the optimal solution of the dual of ( ) MKP then the optimal value of (KP,w*) is a tight upper bound for (QMKP) and it is obtained in a very competitive CPU time.Nevertheless, the direct expansion used in the above upper bound procedure results in an increase of the size of the linearized problem. Indeed, the number of 0-1 variables is equal to ∑ = well known that fewer variables are included in the program so less running time is consumed. The purpose of the next sections is to try to reduce the equivalent linearized problem size, using other expansion techniques and to evaluate the impact of such techniques on the computation of the upper bound

  ,…,u j and j=1,…,n. Since the integer variables are now replaced by 0-1 variables, we transform the resultant 0-1 quadratic program into a 0-1 linear problem as follows: we replace each objective function termf j (x j ) = c j k -d j k² .The problem (QMKP) is thus equivalent to the following problem (MKP 2 ):

=

  For a given index j, let us consider the following linear problem ( ) j P and its dual problem ( )

  continuous knapsack problem: for a given index j, all the variables jk y , coefficients of these variables in the objective function g(y) are in decreasing order (because of the concavity of function last equality holds for the same reason as in the previous paragraph (the coefficients of variables jk y in the unique knapsack constraint are identical and the coefficients of these variables in the objective function are in decreasing order). This implies that:

  the non applicability of both (III) and a linear interpolation to transform (QMKP) into an equivalent 0-1 piecewise linear problem through the use of only

As previously mentioned, the aim of this study is to compare other expansion techniques in the upper bound process developed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Since, the variables are now binary the next step of the algorithm of Quadri et al. [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] Proof 4.1 Assume that it is possible to replace each function f j (x j ) = c j x j -d j x j ², for all j from 1 to n, by a linear function g j (z j ) using exactly ⎡ ⎤ ) ( log 2 j u 0-1 variables. We denote by (H) this assumption.

If (H) is true then it should exist coefficients g jk such that:

The assumption (H) must be satisfied for all problem data. Let us set u j = 3, c j = 10 and d j = 4.

Consequently,

-if x j = 1 (cf. ( 2) and ( 3)) then g j1 = 6 (i), z j1 = 1 and z j2 = 0.

-if x j = 2 (cf. ( 2) and ( 3)) then g j2 = 4 (ii), z j1 = 0 and z j2 = 1.

-if x j = 3 then (cf. ( 3)) then z j1 = z j2 = 1 which implies (cf. ( 2)) that 3c j -9 d j = g j1 + g j2 = -6 (iii).

The equations system (I), (II) and (III) has clearly no solution. Consequently there is a contradiction with (H). Proposition 4.1 shows the non applicability of a binary expansion of the integer variables for (QMKP) so as to transform the initial problem into a 0-1 linear program.

Consequently, the upper bound method proposed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] can not be applied together with expansion (III).

Section 5. Concluding remarks

In this paper we have theoretically compared the use of three techniques to rewrite integer variables into zero-one variables in an upper bound procedure for (QMKP) developed by Quadri et al. [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF], which provides, to the best of our knowledge a bound closer to the optimum than the existing methods. More specifically, we have compared a direct expansion of the integer variables, originally employed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] with a direct expansion with additional constraints (II) and with a binary expansion (III).

We have proved that (II) provides the same upper bound as the one computed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] whereas it involves n additional constraints. We therefore do not expect an improvement of the upper bound computational time. Finally, we provide a proof of