. This work presents the estimation of the parameters of an experimental setup, which is modeled as a system with three degrees of freedom, composed by a shaft , two rotors, and a DC motor, that emulates a drilling process. A Bayesian technique is u sed in the estimation process, to take into account the uncertainties and variabilities intrinsic to the measurement taken, which are modeled as a noise of Gaussian nature. With this procedure i t is expected to check the reliability of the nominal values of the physical parameters of the t est rig. An estimation process assuming that nine parameters of the experimental apparatus ar e unknown is conducted, and the results show that for some quantities the relative deviation w ith respect to the nominal values is very high. This deviation evidentiates a strong deﬁcienc y in the mathematical model used to describe the dynamic behavior of the experimental apparatus . Mecánica Computacional Vol XXXII, págs. 763-773 (artículo completo) Carlos G. García Garino, Aníbal E. Mirasso, Mario A. Storti, Miguel E. Tornello (Eds.) Mendoza, Argentina, 19-22 Noviembre 2013 Copyright © 2013 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar 1 INTRODUCTION A drillstring is a device, used to drill oil wells, whose thre e-dimensional dynamics behavior is very complex, and subject to three mechanisms of vibration: longitudinal, transverse, and torsional ( Spanos et al. ,2003). The torsional mechanism is responsible for causing a type of dynamics called stick-slip. In this mod e, the rotation of the column is blocked due to dry friction with the borehole, stick, and t hen suddenly released, slip, to be afterwards blocked again. So, due to the constant rate of r otation imposed at the top of the column, the column is severely twisted. This may cause fatigue damage, or even the rupture of the structure. To characterize the stick-slip mode and, consequently, to d evelop mechanisms to avoid it, the works of Cayres et al. (2013) andCayres(2013) present numerical studies which investigate, by means of a sensitivity analysis, the inﬂuen ce of the friction models and the caracterization of the structure parameters during the occurrence of the stick-slip. Experimentaldataobtainedfromatestapparatus, inlabora toryscale, areusedtovalidate the mathematical model developed. One of the challenges associated to the use of an experimenta l setup to validate a mathematical model is the correct characterization of the p arameters that are used to describe the behavior of the experimental setup physical sy stem. Commonly, a physical phenomenon is not described perfectly by a mathematical mod el, which may be due to an inability of the model itself, or also, due to inaccuracie s of the measuring instruments, ambient noise, human errors, and several others factors. In this way, the measures that one has for estimate the parameters of a physical system are s ubjected to variabilities ( van der Heijden et al. ,2004). Fortunately, the measurements erros can be circumvented with the use of stochastic techniques of estimation, in part icular, the Bayesian inference, which updates the estimative of the parameters for each new e vidence collected about them (Bolstad,2007). In order to characterize the experimental setup presented i nCayres et al. (2013) and Cayres(2013), that emulates a drilling process, this paper presents a st udy in which the physical parameters associated with the test rig are est imated through the Bayesian technique. The main objective that one wants to achieve with the parametric estimation is to check the reliability of the nominal values of the physi cal parameters of interest. Thispaperisorganizedasfollows. Insection 2theexperimentalapparatusisdescribed. The mathematical modeling of the physical system appears in section3. In section 4, the Bayesian estimation technique is brieﬂy described. Then, i n section 5, the estimation results for the experimental setup are presented and discus sed. Finally, in section 6, the main conclusions are emphasized and some directions for fut ure work outlined. 2 EXPERIMENTAL SETUP Theexperimentalapparatususedintheworksof Cayresetal. (2013)andCayres(2013), which is intended to emulate a drilling process, basically c onsists of a large ﬂexible shaft, one DC motor and two rotors, as can be seen in the schematic rep resentation of the Figure1. In this experimental setup, the ﬂexible shaft has a very low ﬂ exural stiﬀness, so that it emulates the behavior of a drillstring, which is a very sle nder structure due to the high value of the ratio ls/Ds, between the column length, ls, and the diameter of the column cross section, Ds.M.G. SANDOVAL, R. SAMPAIO, A. CUNHA JR 764 Copyright © 2013 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar The DC motor has the role of imposing a constant speed on top of the column, so that it can penetrate the rock. The rotor r2is located between the DC motor and the shaft, and concentrates the inertia of the motor, which is denoted b yjm. On the other extreme of the shaft one can ﬁnd the rotor r1, which has an inertia j1. The modeling of the experimental setup assumes that it has th ree degrees of freedom, being the ﬁrst degree of freedom (DoF) the angle of rotation o f the shaft around the rotor r1, which is denoted by θ1, the second DoF is the angle θ2, which measures the shaft rotation with respect to the rotor r2, and the last DoF is the the electrical charge in the motor, denoted by q. The other physical parameters of this mechanical system, as w ell as their nominal values, are presented in Table 1. For further details about this experimental apparatus and t he model that is used to describe the experiments, the reader is encouraged to see th e work of Cayres(2013). Figure 1: Schematic representation of the experimental setup used t o emulates a drilling process. 3 MATHEMATICAL MODELING The following lines reproduce the mathematical modeling, d eveloped by Cayres(2013) andCayres et al. (2013), of the experimental setup described in the last section. Starting the equationing by the motor, one has that the sum of the electric potentials in the motor is equal to the external potencial applied on it, so that lm¨q(t)+rm˙q(t)+ke˙θ3(t) =v, (1) wherelmis the internal inductance of the motor, rmis the internal resistance of the motor, keis a constant of voltage, and θ3is the angular rotation of the motor axis. Note that the upper dot denotes diﬀerentiation with respect to time. The balance of the torques acting on the motor provides the fo llowing equation jm¨θ3(t)+cm˙θ3(t)−kT˙q(t)+tf=−tm, (2) wherecmis the internal damping of the motor, kTis a constant of torque, tfis the torque due to motor internal friction, and tmis the torque on the motor due to the coupling between the motor and the gearbox.Mecánica Computacional Vol XXXII, págs. 763-773 (2013) 765 Copyright © 2013 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar Column parameters description value unit j1moment of inertia 28 .3×10−3kgm2 lslength of the column 2 .4m Dsdiameter of the column 3 mm mr1mass of the rotor 1 6 .4kg kstorsional stiﬀness of the column 3 .0×10−1Nm/rad rr1radius of the rotor 1 188 mm t1friction torque in the rotor 1 0 Nm Motor parameters description value unit jmmoment of inertia 4 .0×10−4kgm2 lminternal inductance of the motor 1 .1×10−3H rminternal resistance of the motor 3 .3×10−1Ω kTconstant of torque 1 .2×10−1Nm/A keconstant of voltage 601 .6×10−4V/(rad/s) tftorque due to motor internal friction 1 .0×10−1Nm cminternal damping of the motor 18 .0×10−5Nm/(rad/s) imgear reduction factor 1 /8 8:1 vvoltage 8V Table 1: Numerical values for the physical parameters of the experiment al setup. Due to the reduction factor in the gearbox, there is a relatio nship between the angles of rotation θ2andθ3, which is given by imθ3=θ2, (3) beingimthe gear reduction factor. This reduction factor also provi des a coupling between tmand the torque on the second DoF, t2, which is expressed by tm=imt2. (4) Moreover, the balance of the torques which are acting on the r otorr2provides −ks/parenleftbig θ1(t)−θ2(t)/parenrightbig =t2, (5) beingksis the torsional stiﬀness of the shaft.M.G. SANDOVAL, R. SAMPAIO, A. CUNHA JR 766 Copyright © 2013 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar Likewise, the balance of the torques acting on the rotor r1provides the following diﬀerential equation j1¨θ1(t)+ks/parenleftbig θ1(t)−θ2(t)/parenrightbig =−t1. (6) Combining the Eqs.( 2) to (5), one can write jm¨θ2(t)−i2 mks/parenleftbig θ1(t)−θ2(t)/parenrightbig −imkT˙q(t)+cm˙θ2=−imtf. (7) Thus, through a rearrangement of the Eqs.( 1), (6), and (7), it is possible represent the dynamics of the experimental setup by the following state sp ace equation ˙y(t) = [A]y(t)+F(t), (8) where the state vector y(t), the control vector F(t), and the system matrix [ A] are, respectively, deﬁned as y(t) = θ1(t) θ2(t) q(t) ˙θ1(t) ˙θ2(t) ˙q(t) ,F(t) = 0 0 0 −t1/j1 −imtf/jm v/lm , (9) and [A] = 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 −ks/j1ks/j10 0 0 0 i2 mks/j1−i2 mks/jm0 0−cm/jmimkT/jm 0 0 0 0 −inke/lm−rm/lm . (10) Later in this work will be necessary to integrate the initial value problem deﬁned by the ODE of the Eq.( 8), and an appropriate initial condition, assumed as the zero vector for all the simulations reported in this work. This integration is c arried out numerically, using the routine ode23t of MATLAB, that uses a trapezoidal rule with“free”interpolant to construct the time response of the dynamical system under an alysis. The choice of this ODE solver was motivated by the presence of a moderately stiﬀ ness in the system of equations to which one wish to solve. 4 PARAMENTERS IDENTIFICATION The Bayesian statistical inference theory provides the lin k between the results of the experimental observations with the theoretical predictio ns. This can be tapped for para- metric estimation techniques. Thus, given the known distri bution of the observations con- ditioned on the parameters p(y|x) (forward PDF), it is desired to ﬁnd a probability distri- bution of the parameters inﬂuenced by the observations p(x|y) (posterior PDF), ( Bolstad, 2007).Mecánica Computacional Vol XXXII, págs. 763-773 (2013) 767 Copyright © 2013 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar Thus, p(x|y) =1 p(y)p(y|x)p(x), (11) wherexis the parameters vector and yis the sample vector. The Eq.( 11) is the Bayes theorem, and states that the posterior PDF can b e obtained from the forward PDF, determined by the experiment (or simul ation). The a priori PDF, p(x), iscalculatedbasedonwhatisknown(orbelievedtobeknow n)abouttheparameters before observations. The ratio 1 /p(y) is a normalization constant. To obtain the desired parametric estimation through Bayesi an inference, it is necessary to seek the set of parameter that maximize the Eq.( 11). For this purpose, the natural logarithm is applied to the posterior PDF, considering Baye s law and under the assump- tion that the values of the parameters before the observatio n do not change i.e., p(x) is constant. In this way one obtains that ln/parenleftbig p(x|y)/parenrightbig ∝ln/parenleftbig p(y|x)/parenrightbig , (12) that is, the posterior PDF is proportional to the forward PDF . Therefore, to maximize the left hand side of ( 12) is similar to maximize its right hand side. Ingeneral, themeasurementsofanexperimenthavesometype ofdispersion, associated with each observation, and this dispersion can be modeled as an error of random nature, which will be called n. This error is added to the value given by the physical theory ˆy, which is considered free of error. Thus, one has y=ˆy+n. (13) It is assumed that the ﬁrst two moments of nare known (the average value µnand the covariance matrix [Γ n]) and that its distribution is a Gaussian. Also, it is supposed t hat the measurements are independent and identically distribu ted (i.i.d.) events, so that the covariance matrix is [Γ n] =σ2 n[I], where [ I] is the identity matrix and σnis the Gaussian standard deviation. Therefore, the misﬁt function is deﬁned as ε(x;y) =1 σ2 n(y−ˆy)T(y−ˆy) =1 σ2 n∝bardbly−ˆy∝bardbl2, (14) and must be minimized in order to estimate the parameters. Th is procedure is called least-squares method. The procedure followed in order to estimate the parameters o f the experimental appa- ratus consists of collecting the experimental data, constr uct the misﬁt function, Eq.( 14), and then minimize it. The minimum of the function given by the Eq.(14) is searched with the aid of the function fminsearch of MATLAB, which used M-simplex method ( Nelder and Mead ,1965). For further information, the reader can see van der Heijden et al. (2004).M.G. SANDOVAL, R. SAMPAIO, A. CUNHA JR 768 Copyright © 2013 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar 5 RESULTS AND DISCUSSION 5.1 Veriﬁcation of the technique eﬀectiveness First, to verify the eﬀectiveness of the method, only two par ameters are estimated, namely the constant voltage keand the theoretical internal damping of the motor cm. Thus, the parameter vector is written as x=/parenleftigg cm ke/parenrightigg . (15) Theotherparametersvaluesareareassumedtoknown, andare obtainedfromtheTable 1. Given that, in the laboratory, it is only possible to measure the values of the rotations, θ1andθ2, and their respective time derivatives, ˙θ1and˙θ2, a vector of samples is y=/parenleftig θ1θ2˙θ1˙θ2/parenrightigT =ˆy+n, (16) whereˆyis given by the numerical solution of the Eq.( 8), and the noise n, as previously stated, is assumed to be i.i.d. and Gaussian. Thus, the misﬁt function, according to Eq.( 14), is given by ε(x;y) =/parenleftbigg1 σ2 n/parenrightbigg2/summationdisplay i=1/bracketleftigg/parenleftig θi−ˆθi(x)/parenrightig2 +/parenleftbigg ˙θi−ˆ˙θi(x)/parenrightbigg2/bracketrightigg (17) withi= 1,2. In this eﬀectiveness veriﬁcation test, the ˆθi(x) andˆ˙θi(x) are given by the mathematical model, while θiand˙θiare obtained from ˆθi(x) andˆ˙θi(x) by adding a noise. Estimatesofthetwounknownparameters, fordiﬀerentvalue softhestandarddeviation σn, using as initial guess cm= 0.001 and ke= 0.01, and the corresponding deviation, computed with respect to the nominal values of the Table 1can be seen in Table 2. The results reported in Table 2show a good agreement between the values estimated by the Bayesian technique and the reference values shown in T able2. Relative deviations less than 1% can be seen in the case of cm, and negligible values are observed for ke. These results provide an indication that the available nominal va lues forcmandkeare realistic. cm ke σnestimation relative deviation estimation relative deviatio n 0.001 18 .9×10−50.77% 601 .6×10−40.0064% 0.01 19 .1×10−50.77% 601 .6×10−40.0068% 0.1 19 .0×10−50.60% 601 .6×10−40.0051% 1.0 18 .9×10−50.82% 601 .9×10−40.0061% Table 2: Parameters of the experimental setup estimated for diﬀerent values of σn, using as initial guess cm= 0.001 and ke= 0.01.Mecánica Computacional Vol XXXII, págs. 763-773 (2013) 769 Copyright © 2013 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar 5.2 Estimation of the experimental apparatus parameters In this section are presented the results of the estimation o f all the nine physical parameters that are associated to the experimental setup. I n this way, the parameter vector is given as follows x=/parenleftig jmcmkermktlmksj1tf/parenrightigT . (18) For this case the expression of the misﬁt function is similar to the one presented in Eq.(17), but now the sample vectors are no longer calculated based o n the mathematical model, but are obtained experimentally. The ﬁrst thing to notice is that when one tries to estimate all the nine parameters, simultaneously, using the Bayesian approach used in the est imation procedure reported in the last section, convergence problems in the minimizati on procedure are observed. To circumvent these problems of convergence, a heuristic pr ocesudre to minimize the misﬁt function is proposed. In this heuristic algorithm, ﬁr st it is assumed that the pa- rameters x1= (tfjm)Tare unknown, and that the other parameters are given by the values of the Table 1. Then, to obtain estimates for tfandjm, the least-squares method is executed for 10 iterations, using an arbitrary initial gu ess. In the next step of the heuristic algorithm, the unknown vect or is assumed to be x2= (jmcm)T. The value of the parameter tfand the initial guess for jmare given by the estimation of the last step. Again the least-squares method i s executed for 10 iterations, and an estimation for jmandcmis obtained. In what follows, the procedure above continues using the fol lowing order to the vector of unknowns: x3= (cmke)T,x4= (kekt)T,x5= (ktrm)T,x6= (rmlm)T,x7= (lmks)T, x8= (ksj1)T,x9= (j1tf)T. The overall procedure is repeated until a steady-state is re ached for the misﬁt function. A schematic representation of this heuristic can be seen in F igure2. SeeSandoval (2013) for further details. Figure 2: Schematic representation of the proposed heuristic.M.G. SANDOVAL, R. SAMPAIO, A. CUNHA JR 770 Copyright © 2013 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar Despite its heuristic nature, the algorithm shows good conv ergence properties, as can be seen in the Figures 3and4, which respectively show the value of the misﬁt function for each iteration, and the comparison between the rotation s and their time derivatives, obtained from the experimental measurements and from the es timated curves. Note that the good agreement between the estimated and measured curve s increases the degree of conﬁdence in the heuristic estimation. Figure 3: Value of the misﬁt function for each iteration index of the heur istic algorithm. Figure 4: Values measured (blue line) in the experimental setup and curves (red marks) constructed by the Bayesian estimation method.Mecánica Computacional Vol XXXII, págs. 763-773 (2013) 771 Copyright © 2013 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar Estimated values for the physical parameters of interest ar e presented in Table 3, as well as the respective initial guess used in the estimation, and the relative deviations, which are measured with respect to the nominal values shown o f the Table 1. Note that some parameters exhibit a small relative deviation, with a m agnitude smaller than 7%, while others estimates are very distant from the nominal val ue, with a deviation of up to 60%. Note that the nominal values in the Table 1are in fact not the correct value of the parameters of the experimental setup, but, in fact, what one believes to know about them. Besides that, as the heuristic used showed itself reliable, one can argue that the nominal values of the parameters jm,rm, andlmare not close to the real values of the test rig parameters. However, it is the authors belief that this discrepancy in dev iations has another expla- nation. First one can note that all the parameters have large deviations are associated with the motor. Note also that the modeling of the motor assume s a linear behavior, which does not correspond to reality, since torsional dynam ics of the shaft is inherently nonlinear due to its coupling with the motor. Another factor worthy of discussion is the representativenes s of a model with only two degrees of freedom of rotation to described the dynamics of a ﬂexible shaft, which is actually a deformable body with inﬁnite degrees of freedom. Therefore, it is the main conclusion of this study that the de ﬁciencies of the model used to describe the dynamics of the experimental setup are respo nsible for causing the large deviations, with respect to the reference values, observed in the estimated parameters. The authors believe that better estimates can be obtained if a more realistic model is employed. For example, using a system with distributed para meters for modeling the shaft and to include a transistor in the description of the mo tor. parameters initial guess estimation relative deviation jm 4.0×10−45.0×10−428.0% cm 19.0×10−518.0×10−56.6% ke601.6×10−4587.7×10−42.3% kT 1.2×10−11.3×10−15.0% rm 3.3×10−15.3×10−160.0% lm 1.1×10−30.8×10−323.3% ks 2.6×10−12.7×10−13.0% j1 28.3×10−329.8×10−35.4% tf 1.0×10−11.0×10−10.0% Table 3: Parameters of the column-motor physical system estimated by th e Bayesian estimation method. Further results and discussion about the estimation of para meters of the experimental setup of interest can be seen in Sandoval (2013). 6 CONCLUDING REMARKS Inthisworkwaspresentedaprocedureforestimatingthepar ametersofanexperimental apparatus, composed by a shaft, two rotors, and a DC motor, th at emulates a drilling process. A Bayesian approach was used to address the variabi lities that are present in the experimental measurements due to factors such as measureme nt error, noise, etc.M.G. SANDOVAL, R. SAMPAIO, A. CUNHA JR 772 Copyright © 2013 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar For a veriﬁcation test, to check the eﬀectiveness of the Baye sian estimation technique, ﬁrst only two parameters were considered as unknown. The res ults obtained with the estimation procedure for this two parameters problem prese nt a very low deviation from the nominal values of the experimental apparatus for both pa rameters. On the other hand, in a estimation process which considers ni ne parameters of the experimental apparatus as unknown, convergence problems o f the Bayesian method were observed. A heuristic procedure was proposed to overcome th is problem. The estimation results obtained with the proposed algorithm show some para meters with low relative de- viation with respect to the nominal values, and others with l arge deviation. The presence of these large deviations evidentiates a strong deﬁciency i n the model used to describe the dynamic behavior of the experimental apparatus, especi ally by not considering the continuous nature of the shaft and nonlinearities intrinsi c to the shaft/motor coupling. ACKNOWLEDGMENTS The authors are indebted to Brazilian agencies CNPq, CAPES, an d FAPERJ for the ﬁnancial support given to this research. They also acknowle dge Mr. Bruno Cayres and Prof. Hans Weber by providing the experimental data used in th is work.