Wisama Khalil
email: wisama.khalil@ircyn.ec-nantes.fr

Philippe Lemoine

GECARO: A System for the GEometric CAlibration of RObots

Keywords: calibration, geometric parameters, static accuracy, identification, robots. MOTS-CLÉS : étalonnage, paramètres géométriques, précision statique, identification, robots. APII-JESA APII -JESA 2

This paper presents a software package for the simulation and the practical calibration of the geometric parameters of robots. This system which is called GECARO, GEometric CAlibration of RObots, contains a large variety of methods to identify the geometric parameters of robots. GECARO is running on PC computers and developed using MATLAB; any general serial robot can be treated directly. The identifiable parameters are determined using a numerical method based on the QR decomposition, while the identification is carried out using linearized model which is solved iteratively using least squares criterion and by updating the observation matrix after each iteration.

RÉSUMÉ. Cet article présente un logiciel pour la simulation et la réalisation de l'étalonnage de paramètres géométriques des robots. Ce système qui est baptisé GECARO, GEometric

CAlibration of RObots, est composé d'un ensemble de méthodes pour l'identification de paramètres géométriques. GECARO est développé sous MATLAB et s'exécute sur PC ; il peut traiter directement n'importe quel robot à chaîne ouverte simple. Les paramètres identifiables pour chaque méthode sont déterminés par une méthode numérique fondée sur la décomposition QR de la matrice d'observation. L'algorithme d'identification utilise une méthode itérative du type moindres carrés linéaire qui met à jours les paramètres géométriques après chaque itération.

Introduction

The absolute accuracy of a robot depends to a large extent on the accuracy of the values of the geometric parameters used in the direct and inverse geometric models. Geometric calibration of robots is the process by which the parameters defining the base frame parameters, link parameters and end-effector parameters are precisely identified.

Classically, the geometric calibration is carried out by solving a system of linear or non-linear equations which is a function of the joint positions and the location (position and orientation) of the end-effector frame [ROT 87] [WHI 86] [CAE 93] [GUY 95] [DAM 96] [HOL 89]. All the robots are provided with joint position sensors, thus the problem of using the classical methods is to find the appropriate external sensor which can measure with sufficient accuracy the position or location of the terminal frame in the world fixed frame. To overcome this problem, autonomous calibration methods which do not need an external sensor have been proposed [BEN 91] [EDW 94] [ZHO 95] [TAN 94] [KHA 95] [NAK 95], one can find also methods using simple sensors giving the distance between the terminal point of two configurations.

The software package GECARO contains a complete set of these methods, such that we can select the appropriate method as function of the robot and the available sensors.

The paper is organized as follows: the parameters defining the robot will be presented in section 2. In section 3 we describe the calibration methods developed in GECARO. Section 4 presents the software package GECARO. Section 5 presents an example. Section 6 is the conclusion.

Description of the robot

Two types of parameters are needed:

-The geometric parameters defining the location of the end-effector of the robot with respect to the fixed world frame.

-The joint parameters defining the joint limits, the joint offsets, the joint transmission gains, and the coupling matrices between the motor and joint variables.

Definition of the geometric parameters

We consider serial robots consisting of n joints and n+1 links. Link 0 is the base and link n is the terminal link, frame j is defined fixed on link j. We denote:

frame -1: the fixed reference frame, frame n+1: the end-effector frame.

The end-effector location can be calculated with respect to the reference frame by the direct geometric model:

-1 T n+1 = -1 T 0 0 T n (q) n T n+1 [1]
where i T j is the (4x4) transformation matrix defining frame j with respect to frame i, and q is the (nx1) joint variables vector.

The definition of the link frames will be carried out by Khalil and Kleinfinger notation [KHA 86] [KHA 99]. Frame j is defined such that z j is along the axis of joint j, and x j is perpendicular to z j and z j+1 . Frame j is defined with respect to frame j-1 by the matrix j-1 T j , function of the four parameters ( j , d j ,  j , r j ), [Figure 1], such that: j-1 T j = Rot(x, j ) Trans(x,d j ) Rot(z, j ) Trans(z,r j )

[2] Although that frames -1 T 0 and n T n+1 can be arbitrarily defined, it have been proved that the calculation of -1 T n+1 can be obtained as a function of the geometric parameters of n+2 frames represented by four parameters for each frame ( j , d j ,  j , r j ) for j = 0,…,n+1, except for the first frame for which  0 = 0, d 0 = 0 [KHA 91a]

O j-1 j-1 z d x j-1 z j x j O j  j  j r j j
[KHA 99]. The joint variable q j is equal to  j if joint j is rotational and r j if joint j is prismatic. We define the type of the joint by the parameter  j , with:  j = 0 for j rotational,  j = 1 for j prismatic, and  j = 2 if frame j contains only constant parameters as in the case of frame 0 and frame n+1.

If the axis of joint j is parallel to the axis of joint j-1, an additional parameter  j must be considered [HAY 88], the nominal value of  j is equal to zero, such that:

j-1 T j = Rot(y, j ) Rot(x, j ) Trans(x,d j ) Rot(z, j ) Trans(z,r j ) [3]
The calibration object is to identify the deviation of the real parameters from the nominal values, thus to identify:  0 , r 0 , and  j , d j ,  j , r j ,  j (for j = 1, n+1).

Definition of joint parameters

We mean by joint parameters the following:

-Qmax and Qmin: the vectors of joint limits.

-Offset: a vector, containing the nominal values of the offset of joint positions.

-K: the vector of transmission gains [HOL 89]. This parameter can be calibrated in the same manner as the geometric parameters. Thus we can identify also the errors in the transmission gains of the joints K j (for j = 1, n).

-Cm: the (nxn) coupling matrix between the motor variables.

-Cj: the (nxn) coupling matrix between the joint variables.

The joint variables are calculated as a function of motor variables using the following relation:

q = Cj * diag(K) * Cm * qm + Offset [4]
where: q is the (nx1) joint variables vector, qm is the (nx1) vector of motor variables, diag(K) is a (nxn) matrix, whose (j,j) element is equal to K j .

-In addition to the geometric and joint parameters, we need in the calibration process to define a (n+2)x6 matrix which is called the priority matrix. This matrix defines the order of priority of the geometric parameters. Columns 1,...,6 correspond to the parameters , d, , r, , K respectively. The elements of row j give the priority of the parameters of frame j. When calculating the identifiable parameters, the parameters will be ordered as function of the corresponding priority number, the first parameters will be those having the highest priority number. The importance of this matrix comes from the fact that we construct the base of identifiable parameters as those corresponding to the first independent columns. As a general rule we give a high priority to the parameters which can be updated easily in the control system. If the priority of a parameter is equal to zero, it will not be identified.

Calibration methods

The calibration methods of the geometric parameters differ according to the variables used in the calibration model and the measuring devices to carry out them. A unified approach for all of them can be formulated using the framework suggested in [WAM 88] [HOL 96]. According to this formulation the calibration equation can be written in the general form:

0 = f (q,X, r ) [5]
where: X represents the Cartesian variables giving the position and orientation of the terminal frame,  r is the vector of the real (unknown) values of the geometric parameters.

The nonlinear calibration equation can be linearized to get the differential equation:

y = (q,X,)  [6]
where:  =  r -, defines the (Npx1) vector of the errors on the geometric parameters,  is the vector of the nominal values of the geometric parameters.

To identify  r the calibration equation will be applied for a sufficient number of configurations and the global system of equations will be solved using non linear techniques for equation [5] or linear iterative techniques for equation [6].

Calibration using classical methods

This method is based on using a set of configurations for which the joint positions and the corresponding Cartesian coordinates of the terminal frame position or location are given. Beside the joint position sensors, we need an external sensor to measure the position or location of the end-effector.

The location of the end-effector with respect to the reference frame can be calculated using the nominal parameters and the joint positions by the use of the direct geometric model as given in equation [1]. The geometric parameters will be calibrated in order to minimize the difference between the measured end-effector location and the calculated location.

Using a first order Taylor development, a linear differential model defining the deviation of the end effector location due to the differential error in the geometric parameters can be obtained as [KHA 91a]:

x = J(q,)  [7]
where: x represents the (6x1) vector of the position and orientation error (difference between measured and calculated -1 T n+1 ), in Appendix I we show how to calculate x from the elements of the transformation matrices.

J(q,) is the (6xNp) Jacobian matrix of frame (n+1) with respect to the geometric parameters, the calculation of its columns are given in Appendix II.

 defines the (Npx1) vector of the errors of the geometric parameters.

Equation [7] gives 6 linear equations in the unknown vector .

To identify , equation [7] will be applied for a sufficient number of configurations q 1 ,…,q m , the corresponding locations will be measured and x i will be calculated. The resulting linear system of equations will be represented by:

X = W(Q,)  [8]
where:

W =         J(q 1 ) … J(q m ) , X =         x 1 … x m Q = q 1 ,…,q m
W is the observation matrix of dimension (6mxN p ) with 6m >> N p .

It can be seen that if some columns of W are dependent, and W is of rank b, then relation [8] can be reduced to [KHA 91b]:

X = W b (Q,)  b [9]
where W b (q) contains arbitrarily b independent columns of W. The corresponding parameters are defined by the vector  b, they are known as identifiable parameters or base parameters.

Equation [9] will be solved to get the least squares errors solution, given by:

 b = W  b (Q,) X [10]
with W the pseudo inverse of W b .

 b

The geometric parameters will be updated in X and W b . The pseudo inverse solution [10] will be iterated till convergence.

The determination of the identifiable parameters  b must be done before the identification process. They can be obtained numerically using the QR decomposition of a matrix W similar to that defined in equation [8] but obtained using random configurations [KHA 91b]. The outline of this method is given in Appendix III.

The condition number of W b gives an idea about the excitation of the configurations used in the calibration. To get good results with respect to error modelling and measuring noise, the condition number must be close to one [KHA 91a].

If the external sensor gives only the position coordinates of the origin of frame (n+1), then only the first three equations of relation [7] will be used. In this case we must have 3m >> N p . The identifiable parameters must be carried out on the corresponding observation matrix.

Frame link and position link calibration methods

In classical calibration methods, we need an accurate, fast and not expensive external sensor to measure the real end-effector position or location. Many sensors have been proposed in the literature but neither fulfil these three conditions. Frame link and position link methods need only the joint position sensors, therefore they are called autonomous methods.

Position or frame link methods construct the observation matrix using the fact that different configurations of the robot can give the same position or location of the terminal effector [KHA 95]. Thus, if q i and q j represent two different configurations giving the same location of the end-effector frame then:

-1 T n+1 (q i , r ) = -1 T n+1 (q j , r ) [11]
If the values of the geometric parameters in the model  are different than the real values  r , then using a first order development for relation [11] we obtain:

[J(q j , r ) -J(q i , r )]  = x(q i ,q j ,) [12] where: J is the Jacobian matrix of frame (n+1) with respect to the geometric parameters variations, it is defined in equation [7], the expressions of its columns are given in Appendix II.

x is the (6x1) position and orientation differential vector between configurations i and j, that is to say between -1 T n+1 (q i ,) and -1 T n+1 (q j ,), the calculation of x is given in Appendix I. Relation [12] gives 6 equations in the geometric parameters errors.

To estimate the vector  we have to use equation [12] for sufficient number of pairs of configurations to construct a system of equations similar to [8].

If q i and q j give the same position but with different orientations, the observation matrix will be constructed using the first three equations of [12]. This type of equations will be called position link method.

The calculation of the identifiable parameters can be carried out as given in Appendix III, by studying the QR decomposition of an observation matrix W calculated from [12] using sufficient random pairs of configurations q i (k) and q j (k) which give the same end-effector location (or position in the case of position link). Each pair of configurations are obtained by supposing a random configuration q i (k) then calculating the terminal frame location -1 T n+1 (q i ) (or position -1 P n+1 (q i ) in the case of position link) using the direct geometric model, finally q j (k) will be the solution of the inverse geometric model for -1 T n+1 (q i ) (or -1 P n+1 (q i ) in the case of position link) and such that q j  q i . The inverse geometric model is calculated using a general numerical method based on the differential model of the robot.

It is to be noted that the identifiable parameters in these methods are different than those of the classical methods. In fact the following cases can be deduced directly:

-If a column is constant in the Jacobian matrix, then the corresponding parameter will not be identified, because the corresponding column in the observation matrix in this method will be zero. The geometric parameters corresponding to frame 0 and frame 1 verify this case.

-If a column of the Jacobian matrix is equal for all q i and q j , then the corresponding column in the observation matrix will be zero, and the corresponding parameter cannot be identified. For example in the case of frame link method since frame n and frame n+1 are the same for both q i and q j , thus from the expressions of the columns of J, we can deduce that the following parameters are not identifiable: r n ,  n ,  n , d n+1 ,  n+1 ,  n+1 , r n+1 .

The planar calibration methods

In this case the calibration will be carried out using a set of configurations of the robot, whose terminal points are lying in the same plane. This In GECARO two methods are developed, the first makes use of the plane equation, while the second uses the coordinates of the normal to the plane.

The first method: calibration using the plane equation

The general equation of a plane is supposed as:

a x + b y + c z + 1 = 0 [13]
where a, b, c represent the plane coefficients.

As the terminal point (origin of frame n+1) of the robot is in the plane, then:

a Px(q, r ) + b Py(q, r ) + c Pz(q, r ) + 1 = 0 [14] where Px, Py, Pz represent the Cartesian coordinates of the terminal point in the world frame.

In the case of error modelling, the use of a first order development for  and the plane coefficients in equation [14], we obtain:

[Px(q,) Py(q,) Pz(q,) a Jx(q,) + b Jy(q,)

+ c Jz(q,)]       a b c  = -a Px(q,) -b Py(q,) -c Pz(q,) -1 [15]
where: Jx, Jy, and Jz are the first three rows of the Jacobian matrix J defined in [7].

Equation [15] gives a linear equation in , and a, b, c for each configuration.

Sufficient number of configurations must be used to obtain a system of equations similar to [8] which will be solved iteratively to identify , and the plane coefficients.

The coefficients of the plane are initialized by calculating the equation of the nearest plane to the terminal points of the given configurations. The coordinates of these points are calculated using the robot geometric model and the nominal values of the geometric parameters.

The calculation of the identifiable parameters can be carried out numerically, by studying the QR decomposition of a matrix W calculated from [15] using sufficient random points in a given plane.

If the plane coefficients a, b, and c are precisely known, equation [15] can be reduced to contain only  as unknown.

Second method: calibration using normal coordinates

In this method we make use of the fact that the scalar product of the coordinates of the vector normal to the plane and of any vector between two points in the plane

Description of GECARO

GECARO runs on PC computers under Windows 95/98 or Windows NT. A friendly user interface has been developed using C ++ language, while the calculation algorithms have been developed using Matlab 5.2 software. The main page of GECARO is given in [Figure 2], on which we recognize the parameters defining the robot and the following menus: Each of the other menus treats a calibration method. Each of them has the following principal functions:

Identifiable parameters -Points generation -Identification

In the following we detail these functions:

-Identifiable parameters: This function calculates the identifiable geometric parameters (base parameters) for the selected robot and the selected method.

The identifiable parameters can be obtained from the classical parameters by eliminating those which have no effect on the observation matrix or those whose corresponding elements in the priority matrix are equal to zero, and by regrouping the parameters whose corresponding columns in the observation matrix are not independent. These parameters are calculated using QR decomposition of an observation matrix which can be calculated using either random values, or by using real experimental data. This function can also be used to test if some experimental data are sufficiently exciting (giving a good condition number for the observation matrix, and giving the same identifiable parameters as those obtained using random values).

It is to be noted that the identifiable parameters are not uniquely defined, the QR method will give as base parameters those corresponding to the first b independent columns of the matrix W. This function will order the columns of the matrix W and the geometric parameters vector as function of the priority number associated to each parameter. It is more practical to identify, if possible, the parameters which can be updated easily in the control system without changing the symbolic direct and inverse geometric models of the robot. Therefore we use the following rules to define default priorities as input to this function (the user can define arbitrarily his input priority matrix):

-the highest priority is equal to 5 and corresponds to joint offsets and motor transmission gains.

-the parameters r j and d j which are not equal to zero will take the priority 4, -the parameters  j and  j which are not equal to  k/2, with k integer, will take the priority 3, -the parameters of the base location or the end-effector frames, will get the priority 2, -the rest of parameters will get the priority 1.

The user can put equal to zero the elements of the priority matrix corresponding to the parameters which he does not like to identify, either because their values are known accurately, or because he cannot make use of their identified values.

The output of this function is the regrouping relations of the identifiable parameters and a new priority matrix which is similar to the input one, but the priorities of the non identifiable parameters are put equal to zero.

-Points generation: this function generates a file which can be used to test and simulate the functions of the corresponding method. In general the generated file contains two matrices:

Qm: contains the motor variables of the generated configurations, they satisfy the constraint of the selected method.

X: contains the Cartesian coordinates of position or location corresponding to Qm.

The number of configurations is given by the user, this number must be chosen such that sufficient number of equations can be generated.

-Identification: This function gives the solution of the identified geometric parameters and the precision of the obtained solution. Two input files are needed:

-the file containing the priority matrix, -the file containing the input data required to carry out the identification. For instance in the classical method it must contain the matrix Qm of motor variables and the matrix X of corresponding Cartesian coordinates, while in the planar method only a matrix Qm is needed. In the case of position or frame link methods the number of configurations of each set, giving the same Cartesian position and location, will be given in a vector (whose dimension is equal to the number of sets).

Example

In this example we present the identifiable parameters of the Stanford manipulator shown in [Figure 3] and whose geometric parameters are given in [Table 1]. The robot parameters are defined using the menu "Robot" and the different boxes of the main page of GECARO [Figure 2]. The parameters of [Table 1] correspond to 6 joints (frames 1,...,6) and two fixed frames (frame 0 and frame 7). The identifiable parameters of the different methods are given in [Table 2], where: -The parameters denoted by "0" means non identifiable parameters which have no effect on the model, -The parameters denoted by "n" means non identifiable parameters whose effect are regrouped on the other parameters, -The other parameters are identifiable.

1  0 0 1 2 0 2 /   0 2  0.2 0 1 3 1 2 /  0 0 3 r 0 1 4 0 0 0 4  0 0 1 5 0 2 /   0 5  0 0 1 6 0 2 /  0 6  0 0 1 7 2 1.3 0.2 2 /  0.1 0 0 Table 1.
The geometric parameters of the Stanford (RRPRRR) manipulator (units are given in meters for the distances and radians for the angles).

Conclusion

This paper presents the software package GECARO which is devoted for the simulation and the practical calibration of the geometric parameters of robots. The parameters defining the robot, the calibration methods, and the main menus of this package are described in the paper. GECARO is running on PC computers and developed using MATLAB; any general serial robot can be treated directly. This system contains classical calibration methods which require external sensors to measure the terminal link location, as well as autonomous calibration methods which need only joint position sensors. For all the presented methods the identifiable parameters are determined using a numerical method based on the QR decomposition. GECARO can provide for all of developed methods appropriate robot configurations to simulate the capacity and the converge of each method. It has been seen that the identification using a linearized model and by updating the observation matrix and the error vector after each iteration converges in about 5 iterations even if the errors in the geometric parameters are too big (giving up to 60 cm error on the terminal point position).

The future development concerns the calibration of serial robots having nonnegligible elasticity on the joints or on the links and also concerns the identification of the geometric parameters of parallel robots. numerical zero, if the element |R ii |  , the corresponding parameter  i is not identifiable. The numerical zero  can be taken as [LAW 74]:

 = c..max |R ii | [34]
Where  is the machine precision.

To calculate the identifiable (or base) parameters as function of the standard parameters let us permute the columns of W such that the first b columns are independent:

[W.P] = [ ] W 1 W 2 [35]
where: P is a permutation matrix, W 1 represents b independent columns of W, W 2 represents the (c-b) dependent columns of W.

The QR decomposition of (W.P) gives:

[ ] W 1 W 2 = [ ] Q 1 Q 2     R 1 R 2 0 0
Where R 1 is a (bxb) regular matrix. Then it comes:

W 2 = W 1 R 1 -1 R 2 [36]
Thus from [30] and [36] we get:

 = R 1 -1 R 2 [37]

Figure 1 .

 1 Figure 1. Definition of frame j with respect to frame j-1

  idea has been used by different authors [ZHO 95] [TAN 94] [KHA 96] [MAU 96] [IKI 97].

Robot-Figure 2 .

 2 Figure 2. Main page of GECARO

  Figure 3. Stanford manipulator

Acknowledgements

This work has been realized thanks to the financial support of C.R.I.T.T. "Centre Régional d'Innovation et de Transfert de Technologie des Pays de la Loire Productique".

is equal to zero. The main advantage of this procedure is that the coordinates of the normal can be obtained easily using inclinometer [ZHO 95] [MAU 96].

The system general equation is thus: a [Px(q i , r )-Px(q j , r )] + b [Py(q i , r )-Py(q j , r )] + c [Pz(q i , r )-Pz(q j , r )] = 0 [16] where Pu(q i , r ), for u = x, y, z represent the coordinates of the terminal point -1 P n+1 (q i ) at configuration q i . In the case of error modelling, the use of a first order development for  and assuming that the normal coordinates are known, to simplify the writing, we obtain: {a [Jx(q i )-Jx(q j )] + b [Jy(q i )-Jy(q j )] + c [Jz(q i )-Jz(q j )]}  = -a [Px(q i )-Px(q j )] -b [Py(q i )-Py(q j )] -c [Pz(q i )-Pz(q j )]

[17]

where: Ju, Pu are function of q and .

Equation [17] gives a linear equation, for each pair of configurations q i and q j . Sufficient number of configurations must be used to identify . The global system of equations will be solved iteratively and the geometric parameters will be updated after each iteration.

The calculation of the identifiable parameters can be carried out as given in Appendix III, by studying the QR decomposition of a matrix W calculated from [17] using sufficient random configurations where the terminal point of the robot is in a given plane.

Calibration using distance measure

This method needs a set of pairs of configurations where the Cartesian distance between the terminal points of each pair is also given. Thus an external sensor measuring the distance is required. Suppose two configurations q i and q j , where the measured Cartesian distance between the position of their terminal points is D j , thus we can write: r APII -JESA 12 [Px(q j , r ) -Px(q i , r )] 2 + [Py(q j , r ) -Py(q i , r )] 2 + [Pz(q j,  r ) -Pz(q i , r )] 2 = (D ) 2

[18] j r

If the values of the geometric parameters in the model  are different than the real values  r , using a first order development, we obtain:

where: D j is the distance between point j and point i using the nominal model, [Px(q j ,) -Px(q i ,)] 2 + [Py(q j ,) -Py(q i ,)] 2 + [Pz(q j, ) -Pz(q i ,)] 2 = (D ) 2 j Relation [19] can be written under the standard form:

Equation [20] gives a linear equation in . The global system of equations will be solved to get  by using a sufficient number of pairs of configurations.

The calculation of the identifiable parameters can be carried out as given in Appendix III, by studying the QR decomposition of a matrix W calculated using m random points in a given plane.

Calculation of the difference vector between transformation matrices

Let us assume two transformation matrices T i and T j defining frames i and j with respect to the same reference frame such that:

s j n j a j P j 0 0 0 1 , where s, n, a represent the unit vectors along the axes x, y and z respectively.

The vector x representing the position and orientation difference between frames i and j is equal to [KHA 99]:

where: P i -P j is the differential position vector between frames i and j,  i,j is the differential rotation vector between the orientation of frames i and j. It can be calculated by one of the following methods:

-First method:

where x denotes the vector product.

-Second method:

Let the orientation matrix in configuration j be defined as: A j = [ ] s j n j a j , thus:

Where Rot(u,) is the (3x3) orientation matrix representing a rotation  about the vector u.

The differential orientation vector is given by:

Where the values of u and  are obtained from [23] as function of the elements of the matrix A d .

If  is small we can write:

where A(i,j) is the (i,j) element of the matrix A d .

It has been shown that the second method converges better than the first one.

Appendix II

Calculation of the jacobian matrix

Assuming the (4x4) transformation matrix defining frame j with respect to the fixed frame as [KHA 91a] [KHA 99]:

The columns of the matrix J corresponding to the parameters  j , d j ,  j , r j ,  j , k j can be calculated as follows [KHA 99]:

Where:

x denotes the vector product, L j,n+1 is the (3x1) position vector between the origin of frame j and the origin of frame n+1 equal to P n+1 -P j , 0 (3x1) is the (3x1) zero vector, j q j is equal to j  j if joint j rotational or j r j if joint j prismatic.

All the vectors of equations [27] are referred to the measuring fixed frame.

Appendix III

Calculation of the identifiable parameters

Let us consider the following system of overconstrained linear equations:

where W is (rxc) matrix with r>>c. If b is the rank of W, we can write:

where: W 1 represents b independent columns of W, W 2 represents the other (c-b) columns of W.

We can write:

with  is (bx(c-b)) matrix with constant elements.

Using relation [30] , equation [28] can be written as:

with:

The solution of equation [31] will give  b which is called the identifiable (or base) parameters vector. The matrix  is not needed in the identification process.

Numerically, the study of the base parameters is equivalent to study the space spanned by the columns of the matrix W.

Using QR decomposition, W can be decomposed as [DON 79] [LAW 74]:

Q is a (rxr) orthogonal matrix, R is a (cxc) upper triangular matrix, 0 ixj is the (ixj) matrix of zeros.

Theoretically the non identifiable parameters are those whose corresponding elements on the diagonal of the matrix R are equal to zero. Thus assuming  as the


Qussay Jarrar

Farouk Belkadi

Alain Bernard

An Activity-Based Costing Model for Additive Manufacturing

Keywords: Cost model, Cost estimation, Cost drivers, Activity based costing, Additive manufacturing

Additive Manufacturing (AM) has received a high interest in various applications, especially Powder Bed Fusion (PBF) technology. As for any industrial process, cost is one of the most important key performance indicators where good estimation and management have a direct impact on the competitiveness of the enterprise. Therefore, a detailed analysis of the current AM processes cost is essential. This paper has two main contributions, firstly, a critical analysis of the existing cost models in PBF, and metal-based material AM technologies by illustrating their main cost drivers and formulas. Secondly, an Activity-Based Costing (ABC) model is proposed with the aim to cover all important characteristics of AM process. The main cost drivers in AM process are exploited in this model to support the quotation of new product at earlier stages of AM project negotiation. The proposed costing model is part of a global knowledge-based framework for decision aid in AM project.

Introduction

One of the main challenges for decision maker at early stages of AM project is estimating the cost, since it is usually a key point of negotiation [START_REF] Cunningham | Cost modelling and sensitivity analysis of wire and arc additive manufacturing[END_REF][START_REF] Ruffo | Cost estimation for rapid manufacturing-laser sintering production for low to medium volumes[END_REF][START_REF] Busachi | A review of Additive Manufacturing technology and Cost Estimation techniques for the defence sector[END_REF][START_REF] Costabile | Cost models of additive manufacturing: A literature review[END_REF]. Good cost estimation has a direct bearing on the competitiveness of a business of an enterprise because overestimation can result in loss of a business and goodwill market, and underestimation may lead toward financial losses in enterprise [START_REF] Niazi | Product cost estimation: Technique classification and methodology review[END_REF]. Thus, cost estimation Frameworks / models have an important role internally and externally when providing the quotation [START_REF] Siadat | A cost estimation system for manufacturing product using ontology and exepert system[END_REF][START_REF] Mauchand | VCS: value chains simulator, a tool for value analysis of manufacturing enterprise processes (a value-based decision support tool[END_REF]. Estimating the cost is a challenging task which needs a big amount of manufacturing data and knowledge connected to various design and manufacturing aspects: such as material, support, used machine, post-processing, and so on. Face to this complexity, there is recently a high demand of decision support tools and models able to identify the suitable strategy to optimize the global cost while keeping the same level of quality [START_REF] Boenzi | GREENING ACTIVITIES IN WAREHOUSES: AMODEL FOR IDENTIFYING SUSTAINABLE STRATEGIES IN MATERIAL HANDLING[END_REF]. According to the Solid Freeform Fabrication (SFF) symposium held in Austin in 2014 [START_REF]Additive Manufacturing-Developments in training and education[END_REF], most of the experts agreed that future Additive Manufacturing (AM) systems should address functional metallic/multi-material and large-scale parts, allowing low operating cost, speed, and energy consumption. So, understanding the key factors impacting the cost of AM process is necessary for ef-fective deployment and future uptake of high-quality standards [START_REF] Cunningham | Cost modelling and sensitivity analysis of wire and arc additive manufacturing[END_REF]. Indeed, a detailed analysis of the current manufacturing cost and evaluation of the expected improvements reveal a cost reduction potential of about 60% in the next 5 years [START_REF] Facchini | Comparative cost evaluation of material removal process and additive manufacturing in aerospace industry[END_REF], which will significantly advance the market of AM.

In consequence, the costs incurred by AM attract the interest for different stakeholders (i.e., end users, technology developers, AM service providers, investors, etc.). Careful examining of additive manufacturing costs is important to compare additive manufacturing processes with traditional processes, and to identify the main cost drivers used at various steps of the additive manufacturing process [START_REF] Thomas | Costs and cost effectiveness of additive manufacturing[END_REF].

Powder Bed Fusion (PBF) is one of the complex AM processes where thermal energy selectively fuses regions of powder bed to produce parts [START_REF] Standard | Terminology for Additive Manufacturing Technologies[END_REF]. Generally, the process consists of: A thin layer of powder is spread by levelling roller or recoating blade, where a moving energy source (laser beam in case of Laser PBF machines that passes through a system of lenses and reflected by a mirror, used to control the laser beam spot movement on the X, Y planner, onto the platform surface [START_REF] Zhang | Additive manufacturing of metallic materials: a review[END_REF]) melts or sinters the powder into successive cross-sections based on the CAD file. The platform goes down in a specific pitch based on the layer thickness, then the powder is spread for another layer over the solidified layer, to build the next one, until building a complete 3D component. The process may take place in a high vacuum chamber to avoid oxidation issues of metallic powders at high temperature, thus the build chamber is frequently filled with inert gases, argon in most cases [START_REF] Zhang | Additive manufacturing of metallic materials: a review[END_REF]. Examples of PBF process are selective laser sintering (SLS), selective laser melting (SLM), direct metal laser sintering (DMLS), and electron beam melting (EBM) that uses electron beam as an energy source [START_REF] Bahnini | Additive manufacturing technology: the status, applications, and prospects[END_REF].

Despite the significant progress and technology advancement that have been made by PBF, performance in terms of speed, accuracy, process control and cost effectiveness still need to be improved [START_REF] Bhavar | A review on powder bed fusion technology of metal additive manufacturing[END_REF]. This paper has two main contributions. The first one is a critical analysis of the existing cost estimation models in AM, focusing on the ones that are related to metallic materials, and PBF technologies. Then, based on this analysis, the second one is a cost model using Activity-Based Costing (ABC) approach, proposed for PBF laser-based technologies and metal materials, with the ambition to cover some limits of the existing AM cost estimation models. Also, this study represents the main associated cost drivers within the PBF technology process chain. The next section presents a detailed analysis of literature with the aim to identify the main cost drivers and most significant cost models. Then, section 3 sets up the main foundations of the proposed costing framework, while section 4 gives a brief illustration about the application of this framework in an industrial case study.

2

Literature Survey

Cost classification techniques

Cost is the amount that has to be paid to obtain an end product or service regardless of how much the company has gained or lost [START_REF] Huang | The adaptation of product cost estimation techniques to estimate the cost of service[END_REF]. According to the Association for Advancement of Cost Engineering [START_REF]Cost Engineering Terminology[END_REF], cost estimation intends to determine the quantity and predict the costs of constructing a facility, manufacturing goods, or delivering a service. However, there are a wide range of classifications regarding cost estimation [START_REF] Busachi | A review of Additive Manufacturing technology and Cost Estimation techniques for the defence sector[END_REF] based on several aspects like; approach type, granularity level, tools used, and even the application phase. According to Kadir et al. [START_REF] Kadir | Additive manufacturing cost estimation models-A classification review[END_REF], different views of cost models could be categorized into different perspectives depending on the role of the person handling the cost issues like:

Manufacturing perspective: (Task-based classification techniques) includes several phases of product development and manufacturing tasks, it can be grouped into design-oriented (e.g., part design, process planning), or process-oriented that covers cost elements in the production phase, direct and indirect costs (e.g., pre-processing, production (build job), and post processing).

Management perspective: (Level-based techniques) covers a wide range of product costs that are associated with product lifecycle costing, design to cost reduction, remanufacturing and value engineering of the product. These techniques are divided into process-level (costs associated with production, similar to task-based), and system-level (that covers all costs throughout the product life cycle like; maintenance, inventory, etc.).

Financial/accounting perspective: (Method-based) uses a classical classification techniques of cost models based on the used methods and factors. These techniques are categorized into Qualitative (Intuitive techniques, Analogical techniques), and Quantitative (Parametric techniques, Analytical techniques) [START_REF] Niazi | Product cost estimation: Technique classification and methodology review[END_REF].

Intuitive techniques are primarily dependent on past experience (expert knowledge), where this experience can be applied either directly, or using different decision support systems, like (rules, decision tress, case-based, etc.) [START_REF] Niazi | Product cost estimation: Technique classification and methodology review[END_REF]. This approach is suitable at early stages like design, and conceptual stage [START_REF] Busachi | A review of Additive Manufacturing technology and Cost Estimation techniques for the defence sector[END_REF]. Analogical techniques are based on the concept of deriving an estimate from actual information regarding similar real product [START_REF] Ruffo | Cost estimation for rapid manufacturing-laser sintering production for low to medium volumes[END_REF], using historical cost data of products with known cost. Such techniques, like Regression analysis and Neural Networks, have the ability to define a relationship between variables and cost.

The Parametric techniques focus on the characteristics of the product without describing it completely [START_REF] Roy | Cost Engineering: Why, What and How? Decision Engineering Report Series[END_REF] by using the cost estimation relationship that can be presented in mathematical equations, with variables that are associated effectively with cost drivers. The Analytical techniques separate a product into several units, operations, and activities, also the resources consumed during the product life cycle. They express the cost as a summation of all these components [START_REF] Niazi | Product cost estimation: Technique classification and methodology review[END_REF]. It is further classified into different categories (e.g., Operation-based, Feature-based, Tolerance-based, Activity-based, Break-down approach). These are easy and effective methods to apply.

Cost estimation models in AM

The intent of this section is to scrutinize the most relevant cost models defined about the AM "focusing on the ones that are related to metallic materials, and Powder Bed Fusion PBF technologies more specifically". In order to get an overview of the current propositions for cost estimation of AM, in particular the main cost drivers and the used techniques, as well as classify them according to management perspective, described previously, as shown in Table 1.

Process-level cost models

Schroder et al. [START_REF] Schröder | Evaluation of cost structures of additive manufacturing processes using a new business model[END_REF] used Time-driven ABC for the development of a business model which evaluates the process cost of AM technologies (FDM, SLA, SLS, EBM). Thus, their proposition of cost model is consolidated based on detailed subactivities, and an overall of 77 input values are needed for a detailed cost calculation divided into process and order specific information's. Another Time-Driven ABC process-level cost model in AM is proposed by Barclift et al. [START_REF] Barclift | Cost modeling and depreciation for reused powder feedstocks in powder bed fusion additive manufacturing[END_REF]. They provide cost modelling for PBF by expanding the work of [START_REF] Rickenbacher | An integrated cost-model for selective laser melting (SLM)[END_REF], and also defining a new financial depreciation model for reused metal powders using Sum-of-the-years digits depreciation. Their study concludes that the fixed material cost undervalued the cost of build jobs within a range of 13% and 75% for virgin Ti64 powder. The uncertainties in this model are the exact number of reuses permitted for each material alloy, also other cost elements related to quality assurance processes are not included. The work proposed by Cunningham et al. [START_REF] Cunningham | Cost modelling and sensitivity analysis of wire and arc additive manufacturing[END_REF] incorporates the full process chain for Wire-Arc Additive Manufacturing (WAAM) in order to create a detailed cost model. The key cost drivers are identified using sensitivity analysis. The results of this work state that WAAM has a significant potential as cost-effective manufacturing approach compared to other AM technologies. Moreover, Facchini et al. [START_REF] Facchini | Comparative cost evaluation of material removal process and additive manufacturing in aerospace industry[END_REF], provide a cost model to compare the process cost due to production of batch of aerospace parts, adopting both WAAM and traditional machining technologies. Their proposition is using a parametric approach based on computational algorithms, taking into consideration overall manufacturing costs and Non-recurring Engineering NRE costs, which refer to the one-time cost to research, design, develop and test due to production of new parts. Four performance indicators of the process are identified to evaluate the cost of each phase of the process (e.g., setup, building job, part removal, etc.). Ruffo et al. [START_REF] Ruffo | Cost estimation for rapid manufacturing-laser sintering production for low to medium volumes[END_REF] have proposed a cost model based on a parametric approach (as an analytical function identifies the cost estimation relationships) and an engineering approach (to find overall costs by the sum of elementary components used in each step of the production process) in the case of LS (Laser Sintering) manufacturing, as the relationships found are approximations based on statistics. The scope of this cost model is limited to the production process, including material cost as an indirect cost and (overheads, maintenance, capital equipment depreciation, labor, software, and hardware) costs as indirect costs assigned to the components on a machine working-time basis. So, this cost model has to be extended to cover a wider range of AM process chain.

The analogical approach was used in Rudolph et al. [START_REF] Rudolph | Self-learning calculation for selective laser melting[END_REF], and Chan et al. [START_REF] Chan | Data-driven cost estimation for additive manufacturing in cybermanufacturing[END_REF]. Rudolph et al. [START_REF] Rudolph | Self-learning calculation for selective laser melting[END_REF] develop an automated self-learning calculation cost model for SLM implemented within a web-based platform. The customer is able to upload the part geometrical data via online form. Then, the key characteristics (volume, surface and dimensions) are identified to obtain an offer calculation. Linear Regression is used to predict the capacity utilization and the build height in order to estimate build and material costs. The work is only focusing on a cost calculation of the build process including material and manufacturing costs. Chan et al. [START_REF] Chan | Data-driven cost estimation for additive manufacturing in cybermanufacturing[END_REF] propose a cost estimation framework using machine learning algorithms for dynamic clustering, LASSO and Elastic Net Regressions, so the production cost associated with a new job can be estimated based on similar cases in the past. Cost estimation process starts with extracting and calculating 17 different types of features from the process information in Gcode which is generated from the submitted STL file. Those features are then forming the feature vector which will be used as an input of the predictive model. In this study there are some missing links to make effective cost estimation framework, where some cost elements like post-processing, labor cost, and overhead are not included.

System-level cost models

Lindeman et al. [START_REF] Lindemann | Modelling of laser additive manufactured product lifecycle costs[END_REF] aim to provide a deep understanding of AM product lifecycle cost structure for laser-based technologies. The proposed tool is divided into two segments, the first maybe filled by engineer, and the second segment for experienced AM user. It used analytical network process method, multi criteria part classification, then Time-driven activity-based costing as a general method for the calculation, and a good means for allocating cause-based costs to allow the comparison of additive and traditional manufacturing methods. It also, evaluates the relevant factors for specific product, process, warehousing, and capital issues. The cost model is difficult to be implemented in industrial context because of the lack of knowledge about the processes and difficulties for companies to incorporate the lifecycle coasting approach.

Kamps et al. [START_REF] Kamps | Cost-and energyefficient manufacture of gears by laser beam melting[END_REF] carried out qualitative cost analysis incorporating life cycle assessment (LCA) on laser beam melting process. Two integrated models were suggested for cost and life cycle assessment in cradle-to-gate framework focusing on the industrial process sequence. An excel-based cost model is developed for cost calculation as well as build-up duration and argon demand based on company and machine specific data. Whereas the LCA model is based on key assumptions that simplify the model without threaten the calculation outcome significantly. The model needs more experimental measures to improve its accuracy. The results were based for gear production and only few costs were measured. In summary, it is necessary to mention that each of the existing cost models has advantages and disadvantages, but no model meets all criteria satisfactorily, these criteria, among others: simplicity, genericity. For instance, part of the proposed models is limited to AM build process, while others are complex and difficult to implement. So, there is a need for a generic and simple cost model that covers cost elements in different phases of AM process and can be adopted in industrial practices.

3

Proposed cost estimation framework: an ABC approach

ABC method

Activity Based Costing ABC approach has been selected to be used, since it has gained the recognition as an easy and effective method, with better information for controlling capacity cost [START_REF] Huang | The adaptation of product cost estimation techniques to estimate the cost of service[END_REF]. The cost estimation is prepared from the cost incurred by activities during manufacturing of the product [START_REF] Niazi | Product cost estimation: Technique classification and methodology review[END_REF]. ABC could provide the structuring support (the activity) to keep non-financial information such as defect rates (quality), throughput rates (effectiveness of the industrial process) and delivery time [START_REF] Mauchand | VCS: value chains simulator, a tool for value analysis of manufacturing enterprise processes (a value-based decision support tool[END_REF][START_REF] Bosch-Mauchand | Knowledge-based assessment of manufacturing process performance: integration of product lifecycle management and value-chain simulation approaches[END_REF].

It may help the companies to improve product design and manufacturing process, as well as classifying the activities to value-added and non-value added.

The following steps demonstrate the implementation of the ABC approach [START_REF] Ben-Arieh | Activity-based cost management for design and development stage[END_REF]; it starts with identifying the cost centers connected to the directly used resources, including human resources (designer, labour, etc.), and equipment/ machines. Then it analyzes the overall cost (like overhead costs) associated with these cost centers, and it calculates their cost driver's rate. After that it is necessary to assign resources to each cost center, and to determine cost center rates based in the resource cost drivers. Later, the method allows identifying the activities that take place in the product development process, and defining activity driver(s) for each one. Finally, the achievement is finding the total cost of each activity based on the cost-centre resources for each activity, using cost center rate multiplied by the amount of the drivers consumed by each activity (see Fig. 1).

Main Cost Drivers

The activities of the PBF production process chain were identified according to the major production processes involved that represent the scope of this cost model. One of the main purposes of cost modeling is to identify the factors affecting the cost during various processes of AM which are called cost drivers. Thus, before starting with the equations of the cost model, the main cost drivers within the process chain -that are needed in order to estimate the product cost -are illustrated in Fig. 2. It is also important to consider the various interdependent aspects like; product characteristics, customer specifications, machine/ material types, and manufacturing process parameters, in terms of their correlation with the main cost drivers therefore to increase the robustness of cost investigation.

Build time is one of the most essential factors to estimate the cost of AM. Certain machine/ material combination have strong impact on build time. Also, part characteristics, and requirements (part volume, part features, part quality), and manufacturing parameters (hatching and contour percentage, layer thickness, scanning speed, build height, support required) affect the build time. The estimation of build time is beyond the scope of this paper. However, several software packages exist for estimating the build time like (Materialise-Magics, EOSprint, AddUp Manager, etc.). Besides, several works have been caried out for estimating the build time using Machine learning techniques like [START_REF] Munguía | Neural-network-based model for build-time estimation in selective laser sintering[END_REF][START_REF] Angelo | A neural network-based build time estimator for layer manufactured objects[END_REF][START_REF] Mccomb | Predicting Part Mass, Required SupportMaterial, and Build Time via Autoencoded VoxelPatterns[END_REF].

Fig. 2. Main Cost Drivers in PBF process chain

Machine/ equipment cost is also to be taken into account, especially AM machine cost due to the high purchase price of AM machines. This cost is indicated as hourly rate for running the machine (cost/hour). It is based on the purchase price of AM system, service year, and other maintenance expenses associated with such system. Specific machine component like build platform has a separate cost because in some cases, it might have a significant effect on the cost of AM products, this type of cost is based on the dimension "mainly the thickness" of the platform, platform material type, and platform state if it is new or has been already used before for several builds.

Material cost has to be considered carefully because the cost of AM metal material is significantly high compared to raw material used for traditional manufacturing. Two things have to be considered in determining this cost. On one hand, the material price per kilogram which is based on the material type, material supplier, and material state (new or reused). On the other hand, the quantity of the consumed material which is related to the total volume of part, other test parts "that might be included in the build job", support structure, number of parts, and the material waste that includes both material loss during the build process and trapped material within build parts.

Labor cost is related to the labor needed to execute a task or activity (designer, AM operator, post-processing operator, etc.) such as job preparations, machine output, post-processing, and control process. This cost is indicated as an hourly rate cost.

Post-processing and quality control costs may take up to 50% of the final product cost [START_REF] Lindemann | Modelling of laser additive manufactured product lifecycle costs[END_REF]. They both depend on the customer requirements, part characteristics, customer activity sector (function of the product), required labor, used technique, used machine/ equipment, and the duration/ time needed for the activity, which are the most related factors for these costs. Post-processing costs are the sum up of all necessary post-processing operations after AM machine output, (e.g., powder removal, separation of the part from the platform, support removal, machining process, and heat treatment, etc.). On the other hand, quality control cost is the cost of different steps needed to assure and evaluate the various aspects to achieve part quality, and product requirements (e.g., 3D scan, CMM, Tomography, Microstructural evaluation, Metallurgical test, etc.). Indeed, it is substantial to mention that these processes have to be considered beforehand the manufacturing. In specific situations, some of these processes are carried out by subcontractors.

Overhead costs are assigned to the operations of AM system and categorized into: Production overheads, like building rent, support equipment, or energy costs; Administrative overheads that are related to computer equipment, communication, and software license costs. Some cost models (e.g., [START_REF] Lindemann | Modelling of laser additive manufactured product lifecycle costs[END_REF][START_REF] Ruffo | Cost estimation for rapid manufacturing-laser sintering production for low to medium volumes[END_REF]), have included overhead costs as indirect costs along the machine cost, by taking the annualized overhead costs over the machine working time per year. Ill-structured costs are the hidden costs in the supply chain [START_REF] Thomas | Costs and cost effectiveness of additive manufacturing[END_REF] such as, inventory, transportation, supply chain management (purchasing, operation, distribution, integration), and quality (in terms of process failure during machine operation and part rejection) costs. These costs are difficult to consider and not well understood due to the limited knowledge and data in the accounting practices [START_REF]Additive Manufacturing-Developments in training and education[END_REF]. However, AM has a positive potentiality in reducing these costs compared to conventional manufacturing technologies, because of its ability to manufacture parts on demand, produce an entire product in the same build at once, and reduce the link in supply chain by bringing the manufacturers closer to consumers.

Cost calculation for AM

In this section part of the proposed cost model equations are presented. The estimation of AM final product cost is carried out by identifying at first the main needed activities to realize the product, then the cost centers associated with these activities and their drivers. In most of the activities the main drivers are the duration/ time re-quired to achieve the activity. So, this cost estimation model could be also included as a Time-Driven ABC.

CAD/CAM preparation cost = T CAD/CAD prep * (C Designer + C WS)

T CAD/CAD prep: the time required for creating support, simulation, etc.

C Build machine = T Build * C AM-Machine (3)
C Build machine: the cost for operating AM machine (€) C production-stop: the cost of stopping the production (not failures), but to fill the dispenser, empty the re-coater (mainly related to manual tasks) (€) T Build: the time for building the entire job for specific project (h) C AM-Machine: AM machine hourly rate (€/h)

C Build material = M material * C material (5)
M material = (V Build * P melted-density) + M waste [START_REF] Siadat | A cost estimation system for manufacturing product using ontology and exepert system[END_REF] C Build material: the cost of used powder in AM builds job process (€) M material: mass of used material in build job (kg) C material: the commercial price of the used powder (€/kg) V Build: total build volume for specific project (mm3) P melted-density: powder melted density (around "Powder density/0.6") (kg/mm3) M waste: powder loss during AM builds job (kg)

Machine output cost = T rem * C op [START_REF] Mauchand | VCS: value chains simulator, a tool for value analysis of manufacturing enterprise processes (a value-based decision support tool[END_REF] T rem: the total time required to remove the entire build, and cleaning (h) Other costs related to post-process activities are also included in the cost estimation model like; heat treatment, support removal, surface finishing, and machining. They are mainly based on the time/ duration to realize the activity, post-process operator hourly rate, and used machine/ equipment cost rate. Different quality control activities are included as well, parts of these equations are presented below in (11,12, and 13).

Powder removal cost = T PR * (C op + C PR machine) (8)
Project Monitoring cost = T follow up * C project leader [START_REF] Thomas | Costs and cost effectiveness of additive manufacturing[END_REF] T

T QC: time required to realize the control check (h) C QC tool: hourly rate for the used tool/ or machine (e.g., 3D scan, CMM) (€/h) C QC op: quality control operator hourly rate (€/h) Thus, the final equation to estimate the cost of a product is shown below;

Final Product cost = Total material cost + B. Platform cost + ∑Activity cost [START_REF] Zhang | Additive manufacturing of metallic materials: a review[END_REF] Total material cost = C Build material + C PR (Material trapped) + C Psq [START_REF] Bahnini | Additive manufacturing technology: the status, applications, and prospects[END_REF] Total material cost: the total material required for a specific project (€) B. Platform cost: the cost for using a specific build platform (€) C Psq: the cost of the powder sample needed for powder test analysis, if required (€) ∑Activity cost: the summation of the required activities costs (€)

Case Study

To illustrate the application of the proposed cost model, one part designed and produced by laser-based PBF has been selected from [START_REF] Barclift | Cost modeling and depreciation for reused powder feedstocks in powder bed fusion additive manufacturing[END_REF], see Fig. 3. This part was built using EOS Titanium Ti64 material, and manufactured by EOS M280 DMLS machine. The cost will be estimated following the carried-out process presented in [START_REF] Barclift | Cost modeling and depreciation for reused powder feedstocks in powder bed fusion additive manufacturing[END_REF]; the process starts with CAD/CAM preparation by importing the STL file into Materialise Magics software, and Machine setup for the build job, then AM process. Once the build job is completed the entire build is removed. Next, heat treatment and wire EDM (subcontracted processes) are applied. Finally, the required post-processing in terms of support removal and other mechanical process are carried out by the operator. Using the cost model constants from [START_REF] Barclift | Cost modeling and depreciation for reused powder feedstocks in powder bed fusion additive manufacturing[END_REF], the costs for producing the part are estimated using the proposed cost model equations, and compared with the ones estimated in [START_REF] Barclift | Cost modeling and depreciation for reused powder feedstocks in powder bed fusion additive manufacturing[END_REF], as presented in Table 2. The costs are different in both machine setup, and machine output activities because in the proposed model the machine hourly rate cost is not included in these activities, only operator cost rate is included. Material cost includes the fused material to build the part (based on build volume), as well as the powder loss during the build process and the trapped powder inside the build (we assume the same quantity of powder loss from [START_REF] Barclift | Cost modeling and depreciation for reused powder feedstocks in powder bed fusion additive manufacturing[END_REF]).

Discussion and limitation

In comparison with the current propositions for cost estimation in AM (as shown in Table 1), this model is considered as a process level cost model, using ABC approach, in which it has been developed based on an exhaustive review of these propositions, as well as several discussions with the industrial partners. This is done in order to overcome the missing cost elements in previous cost models, and to present the main AM cost drivers as well as their correlation with customer specifications and product characteristics.

The cost model equations cover a wider range of cost elements within AM process chain. In the previous models, the costs for following and creating the manufacturing report for a specific project were not mentioned. In reality, these costs might be significant for complex projects that required, for example, a number of subcontracted activities. Other costs were also highlighted in this paper which are related to control check process (e.g., 3D scan), and to the needed test analysis such as: powder test, metallurgical test, tensile test, etc. The cost for the used build platform was also included since it can vary based on its dimensions and material type. For instance, in some project the customer asks for a specific build platform. Moreover, in this model the quantity of the consumed material was based on several aspects along the part volume such as: the associated test parts, powder loss within the build process, powder trapped inside the part channels, and the powder quantity used of test analysis in some cases.

The cost estimation model could serve as a basic element in two important applications. First it can be adopted in real industrial practices and linked to expert knowledge using different approaches (e.g., rule-based) to support the estimation of cost at early stages. Second, it can be used to track the incurred cost of each activity and its related factors toward a full cost traceability model.

On the contrary, this cost model did not consider the state of the used powder (if it is new, recycled, or mixed). Therefore, the material price is treated as an input for the model, based on the commercial selling price from the supplier regardless the state of the powder. In addition, the model is limited to PBF metal laser-based application. However, the use of this cost model in different applications could be achieved after some updates and modifications. Finally, Costs for build failure and re-design or engineering of pre-existing part were not captured in this model.

Conclusion and Perspective

In this work several cost models in the field of AM have been reviewed, and a cost model using Activity-Based Costing approach is developed. At first, the main Cost Drivers within the PBF process chain have been presented and explained. After that, generic cost estimation equations using ABC are proposed that cover several cost elements like post process, quality control and assurance processes, test analysis process, powder loss, build platform, etc. A case study was presented using Ti64 material part produced by EOS M280 DMLS, in order to clarify the employment of cost model equations. The cost model has been tested with one of our industrial partners in order to estimate the cost based on customer specifications, using developed Business Rules (know-how knowledge) that link customer specifications and product characteristics to cost model inputs and main cost drivers. The proposed cost model is a key step on an ongoing work concerning the development of a knowledge-based decision aid system dedicated to support the whole value chain associated to additive manufacturing project. Several future improvements on this model could be carried out, by further investigation of overhead and Il-structured costs (such as those associated with build failure, build rent, energy consumption, and administrative overhead), on the other hand extending the model by considering other associated costs with product lifecycle from product concept until the disposal of the product.

Fig. 1 .

 1 Fig.1. Activity-Based Costing approach implementation, adapted from[START_REF] Ben-Arieh | Activity-based cost management for design and development stage[END_REF]

  (h) C Designer: The Designer hourly rate (€/h) CWS: The cost rate of using workstation (including software and computer) €/h Machine set-up cost = (T setup + T mat.change) * C op (2) T setup: the time required to prepare and setup the build job (h) T mat.change: the time required to change the material and clean dispenser (h) C op: AM operator hourly rate (€/h) Build job cost = C Build machine + C production-stop

T

  PR: the time required for automatic powder removing (h) C PR machine: hourly cost rate for powder removing machine or equipment (€/h) C PR (Material trapped) = M material PR * C material (9) C PR (Material trapped): the cost of powder lost during the powder removal process, which is the trapped powder within the build parts (€) M material PR: mass of lost powder during powder removing process (kg) Separate parts from platform = (C Sep mach * T Sep mach) + (C Sep op * T Sep op) (10) C Sep mach: hourly rate cost for the used machine (e.g., EDM, Saw) (€/h) T Sep mach: the time/ duration for the separation machine (h) C Sep op: hourly rate cost for the separation process operator (€/h) T Sep op: the time/ duration for the operator (input, remove, monitor) (h)

Fig. 3 .

 3 Fig.3. Automotive Upright, source[START_REF] Barclift | Cost modeling and depreciation for reused powder feedstocks in powder bed fusion additive manufacturing[END_REF]

Table 1 .

 1 Comparison of main AM cost models

	Cost model	Main cost Drivers	Sys-	Pro-	Used meth-
			tem-	cess-	od/approach
			level	level
	Schroder et	CAD preparation cost, operating		✓	ABC (Time-
	al. [20]	cost (preparation process), mate-			Driven-
		rial cost, machine cost, labour			ABC)
		cost (production), post-processing
		cost, quality control process cost,
		administrative cost
	Rudolph et	CAD preparation cost, set-up		✓	Analytical,
	al. [23]	cost, operating cost, labour cost			and analogi-
		(production), production over-			cal "linear
		head, pricing mode cost, material			regression"
		cost, part removal cost, treatment

  follow up: the total time needed to trace and monitor specific project (h) C project leader: project leader hourly rate

	Creating Manufacturing Report cost = T M.rep * C project leader	(12)
	T M.rep: total time required to create the manufacturing report (based on the level of
	this report, and customer specifications) (h)
	Control Check = T QC * (C QC tool + C QC op)

Table 2 .

 2 Compare the estimated costs

	Cost Elements	Estimated costs ($)	Costs from [21] ($)
	CAD/CAM Preparation	630	630
	cost
	Machine Setup cost	550	850
	Build job (C Build ma-	3850	3850
	chine) cost
	Machine output	330	510
	Heat Treatment cost	Stress-relief "Subcon-	350
		tracted" 350
	Separate parts from plat-	Wire EDM "Subcontract-	200
	form cost	ed" 200
	Post Process cost	480	480
	Material cost	2652	2652

Acknowledgement. The presented results were conducted within the French national project "SOFIA" (SOlution pour la Fabrication Industrielle Additive métallique). This project has received the support from the French Public Investment Bank (Bpifrance) and the French National Center for Scientific Research (CNRS). The authors would like to thank all industrial and academic partners for their involvement in this research.