Asma Chargui

Raya El Beainou

Alexis Mosset

Joseph Gavoille

Pascal Vairac

Sébastien Euphrasie

Nicolas Martin

Anisotropic thermal conductivity of nanocolumnar W thin films

Keywords: Glancing Angle Deposition, thermal conductivity, anisotropy, tungsten film

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The GLancing Angle Deposition (GLAD) technique has turned out into an attractive strategy to produce original surface morphologies, especially for designing anisotropic architectures [START_REF] Robbie | Advanced techniques for glancing angle deposition[END_REF]. GLAD films exhibit original physical properties (electrical, thermal, mechanical, and optical behaviors), which are different from the inherent characteristics of the corresponding bulk materials creating new effects [START_REF] Hawkeye | Glancing Angle Deposition of Thin Films[END_REF]. As one of the most interesting physical properties of GLAD films, the thermal conductivity has, until now, been rarely investigated [START_REF] Plawsky | Engineered nanoporous and nanostructured films[END_REF]. This gap of knowledge is connected to the difficulty to perform accurate thermal measurements in thin films, which becomes an even more challenging task in porous and nanostructured media like GLAD films. In addition, the rare studies 2 reporting thermal investigations of GLAD films exhibit thermal conductivity lower than a few Wm -1 K -1 for oxide materials mainly [START_REF] Hass | Low thermal conductivity vapor deposited zirconia microstructures[END_REF]. To the best of our knowledge, no studies have been focused on thermal and electrical properties of GLAD metallic thin films.

In this letter, we investigate thermal conductivity and anisotropy of GLAD W films sputter-deposited with 2 deposition angles: α = 0° and 80°. The 3ω experimental technique [START_REF] Cahill | Thermal conductivity measurement from 30 to 750 K: The 3ω method[END_REF] is used to determine the thermal conductivity of the films, especially the components of the tensor following the two perpendicular in-plane conductivities, namely kx and ky. Surface and cross-section microstructures are observed to correlate the porous character and the anisotropic architecture of the W GLAD films with their very low thermal transport behavior. Deviations from the Wiedemann-Franz law are also discussed considering similarly, the in-plane electrical conductivities (σx and σy.) of the films in connection with their inhomogeneous morphological features.

Material and Methods

W films were deposited on (100) Si and glass substrates by DC magnetron sputtering in a 40 L homemade chamber at a base pressure below 10 -5 Pa. More details can be found in the Supplementary Information. Briefly, an argon sputtering pressure of 0.4 Pa was used for all depositions. No external heating was applied during the growth stage, and depositions were carried out at room temperature. The GLAD technique [START_REF] Cahill | Thermal conductivity of a-Si:H thin films[END_REF] was implemented to produce normal (α = 0°) and tilted columnar architectures (α = 80°). The deposition time was adjusted to get a constant film thickness of 350 nm. The films morphology (top and cross-section views) was observed by scanning electron microscopy (SEM) with a JEOL JSM 7800 field emission microscope.

The thermal conductivity of the films was measured by the 3ω method [START_REF] Cahill | Thermal conductivity measurement from 30 to 750 K: The 3ω method[END_REF][START_REF] Cahill | Thermal conductivity of a-Si:H thin films[END_REF][START_REF] Bogner | Differential 3ω method for measuring thermal conductivity of AlN and Si3N4 thin films[END_REF] (to see the Supplementary Information for experimental conditions). A SiO2 layer (250 -450 nm thick; acting as an electrical insulator and reducing the air/film roughness) was first deposited on W films by ICPECVD (inductively-coupled plasma chemical vapor deposition). A metallic Pt strip (100 nm thick) was then evaporated and patterned by photolithography and lift-of techniques, acting as heater and sensor. The length of the strips was 3 mm for α = 0° and 1 mm for α = 80°, and the width ranging between 2 µm and 30 µm depending on the sample. Strips were either aligned following the x axis (parallel to the particle flux), or the y axis (perpendicular to the particle flux). Isotropic thermal conductivity was calculated for films deposited on Si (α = 0°) using the classical method of the up-shift of the in-phase ΔT line according to the logarithm of the frequency [START_REF] Cahill | Thermal conductivity of a-Si:H thin films[END_REF]. A sample with the SiO2 layer but without the W film was assumed as a reference.

For samples on glass (α = 80°), the model given by Borca-Tasciuc et al. [START_REF] Borca-Tasciuc | Data reduction in 3ω method for thin-film thermal conductivity determination[END_REF] based on a two-dimensional heat conduction model was used to extract in-plane thermal conductivities of the films (kx or ky depending on the strip orientation). For stripes with the smaller width (2 µm), the contribution of the in-plane conductivities becomes important, i.e., ky for strips aligned to the x direction (or kx for y direction). It is also worth noticing that a substrate with a small thermal conductivity (e.g. glass) contributes to the in-plane conductivities. For each sample, the thermal conductivity tensor was extracted from the best fit of doublet for all measurements (different sizes of the strip). The cross-plane thermal conductivity (kz) was calculated but without meaning and thus not given since it includes unknown values of interface resistances.

Electrical conductivity measurements of W films deposited on glass substrate were performed in air at room temperature. Based-on the van der Pauw configuration (4 probe method), Bierwagen et al. [START_REF] Bierwagen | Mobility and carrier density in materials with anisotropic conductivity revealed by van der Pauw measurements[END_REF] proposed an original procedure to determine the in-plane components of the electrical conductivity tensor (and its anisotropy defined as the ratio between electrical conductivity measured following parallel and perpendicular directions to the sputtered particle flux, i.e., σx and σy, respectively) in the perpendicular directions (x, y) of the films.

Results and Discussion

Top and cross-section views by SEM of W films clearly show contrasted morphologies between the two deposition angles (Fig. 1). Conventional sputtered films (α = 0°) exhibit a dense microstructure (top view in Fig. 1a) with homogeneous columns perpendicular to the substrate surface (film's thickness and silica layer on top of the film can be seen in Fig. 1c). W films prepared by GLAD (α = 80°) display an inclined columnar and nanostructured structure (Fig. 1d). Such sputtering conditions lead to the formation of elongated shaped columns following the direction perpendicular to the incident flux (y direction in Fig. 1b). This morphological asymmetry correlates with the shadowing mechanism and a favored diffusion of the incoming sputtered particles in the direction perpendicular to the incident flux [START_REF] Chargui | Influence of thickness and sputtering pressure on electrical resistivity and elastic wave propagation in oriented columnar tungsten thin films[END_REF]. The column tilt angle β, defined as the angle between the substrate normal and the column center axis, reaches 39±2° for these sputtering conditions. It is worth noting that the deposition of SiO2 efficiently covers the surface porous microstructure (surface roughness at the film/air interface is reduced) and also ensures an insulating effect.

Thermal conductivity measurements are presented in Table 1. W films sputter-deposited by the conventional sputtering process (α = 0°) show no in-plane anisotropy. They exhibit a thermal conductivity of 71±5 Wm -1 K -1 at room temperature, which is 41% of the W bulk value (174 Wm -1 K -1 ). Such a difference between thin films and bulk materials is systematically observed. It is connected to the deposition conditions, the film growth techniques and mainly assigned to structural differences. Micro-to nano-crystalline structures reduce thermal conductivity due to an increased phonon scattering from lattice imperfections, point defects or grain boundaries. In addition, the film/substrate interface has to be considered as a barrier to phonon transport, especially in amorphous or poorly crystallized films.

A substantial drop of thermal conductivity is obtained for W films prepared with α = 80° (column angle β = 39°) with kx = 3.2 Wm -1 K -1 and ky = 4.8 Wm -1 K -1 (Table 1). This strong reduction of thermal conductivity is related to the voided architecture, which becomes relevant when the deposition angle α becomes grazing (α > 70°). For α = 80°, it is commonly reported that the films porosity ratio can be higher than 60% [START_REF] Ott | Morphological characterization and porosity profiles of tantalum glancing-angle-deposited thin films[END_REF], which corresponds to a voided structure and thus, a poor thermal conductive medium, i.e., in the same order of magnitude of polymers or even lower.

Imperfections and growing defects are often encountered in GLAD films [START_REF] Beainou | A 4-view imaging to reveal microstructural differences in obliquely sputter-deposited tungsten films[END_REF]. Such defects also reduce the thermal conductivity by lowering the effective mean free path of electrons and phonons. In addition, and due to the high porous structure, films are also partially oxidized, lowering their thermal conductivity. It is important to note that interfacial resistances are not assumed in our model, and thus the out-of-plane conductivity kz is not investigated.

Electrical conduction is similarly influenced by the deposition angle (Table 1). As previously noticed for thermal properties, electrical conductivity of conventional W films (α = 0°) is lower than the bulk value due to electron scattering by surfaces and grain boundaries. For both x and y directions, σ300K is close to 7×10 6 Sm -1 (Table 1). Films prepared by GLAD become more resistive and for α = 80°, σ300K decreases to 1.7×10 5 Sm -1 for y direction. It is even more reduced following x direction with σ300K = 8.4×10 4 Sm -1 . This drop is mainly assigned to the porosity especially favored by the GLAD process, decreasing the electron mean free path. This porous architecture is particularly marked is the direction of the incident particle flux (x direction) and well correlates with the lowest electrical conductivity.

Results in table 1 likewise illustrate the in-plane anisotropy defined as the thermal conductivities ratio (ky/kx = 1.5 for GLAD W films). The columns cross-section of GLAD W films exhibits an elliptical shape with a major axis following the y direction, which favors the heat flow. Columns tend to chain to each other in this y direction whereas the number of pathways available following the x direction is reduced due to a more voided medium. Such anisotropic behaviors have similarly been observed for other physical properties like electrical conductivity and surface acoustic waves [START_REF] Beainou | Electrical conductivity and elastic wave propagation anisotropy in glancing angle deposited tungsten and gold films[END_REF]. These directional-dependent characteristics have systematically been attributed to the fan-like columnar features broadening in the direction orthogonal to the particle flux.

It is worth connecting thermal to electrical conductivity via the Wiedemann-Franz (WF) law, as commonly checked for bulk metals. W films deposited with α = 0° support this law with a Lorenz number L300K of 3.4×10 -8 WΩK -2 whatever the x or y direction. This value completely agrees with the W bulk value (Table 1). Such a WF law is not well satisfied for GLAD W films since L300K = 1.3×10 -7 and 9.4×10 -8 WΩK -2 for x and y directions, respectively. This discrepancy of Lorenz numbers between GLAD and conventional films (and thus bulk) can be assigned to the different effects coming from the free electrons scattering mechanisms (and also to the deposition of SiO2 layer favoring the W oxidation, especially for GLAD films). Transport properties of electrons serving as charge carriers and heat carriers are scattered by grain boundaries and structural defects, particularly present in GLAD films. Assuming the same probability for electrons to be reflected by such grain boundaries and defects, it leads to different effects on thermal and electrical conductivities and thus, to the violation of the Wiedemann-Franz law [14]. It is also worth noting that the highest Lorenz number is obtained following the x direction. This is directly connected to the elongated pores produced in the direction of the column growth (Fig. 1b). Considering that the thermal conductivity strongly depends on the film morphology and volume fraction of pores, their orientation (orthogonal to the particle flux) mainly impedes heat flow according to the x direction, which leads to a low thermal conductivity.

Conclusion

In-plane thermal conductivity of W thin films 350 nm thick was investigated by the 3ω method. Films were sputterdeposited by the GLAD technique using two different deposition angles: α = 0° and 80°. Films with α = 0° showed the expected reduction of the thermal conductivity compared to the bulk, due to the polycrystalline structure of deposited materials. For W GLAD films prepared with α = 80°, the thermal conductivity was two orders of magnitude lower than the bulk one. This strong difference was mainly assigned to the voided architecture typically obtained in films prepared by oblique angle deposition at grazing angles (α > 70°) and their partial oxidation. A correlation with the electrical conductivity well agrees with the Wiedemann-Franz law for conventional W films (α = 0°), whereas a Electrical conductivity (Sm -1 ) 1.9×10 7 (6.9±0.1)×10 6 (8.4±0.1)×10 4 (1.7±0.1)×10 5 Thermal conductivity (Wm -1 K -1 ) 174 71±5 3.2±0.5 4.8±0.5

Lorenz number L300K (WΩK -2 ) 3.1×10 -8 (3.4±0.3)×10 -8 (1.3±0.2)×10 -7 (9.4±1.5)×10 -8

[

  Fig. 1. Top (without SiO2 layer) and cross-section (with SiO2 layer) views by SEM of 350 nm thick W thin films sputterdeposited by DC magnetron sputtering: a) and c) with a normal incidence of the particle flux (α = 0°); b) and d) with a deposition angle α = 80°. The SiO2 insulating layer is also indicated and the blue arrows show the incoming particle flux. The scale bar is the same for all pictures.

Table 1 :

 1 Thermal and electrical conductivities following x and y directions, and corresponding Lorenz numbers at 300 K of W thin films 350 nm thick sputter-deposited with deposition angles α = 0° and 80°.

Corresponding author: nicolas.martin@femto-st.fr

Acknowledgements

This work was supported by the EIPHI Graduate School (contract ANR-17-EURE-0002), the Region of Bourgogne Franche-Comté and the French RENATECH network.

significant difference was clearly brought to the fore for GLAD W films. This discrepancy was attributed to electrons scattering at grain boundaries and by defects, producing different effects on thermal and electrical behaviors. In addition, oxidation favored by the porous structure and the SiO2 layer, also contributed to this divergence since W GLAD films becomes reactive towards oxygen and thus deviate from pure metals.

The thermal anisotropy of GLAD W films was finally determined to focus on the in-plane conduction effect. The inplane anisotropy reached 1.5 due to the elliptical shape of the columns cross-section in the direction perpendicular to the particle flux. This elongated feature of the columns cross-section helps to the heat flow propagation and corroborates other anisotropic behaviors previously reported for electrical conduction and elastic wave propagation.

CRediT authorship contribution statement

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Nahashon O Osinde
email: nahashon.osinde@femto-st.fr

Maxime Etiévant
email: maxime.etievant@femto-st.fr

Jean Bosco Byiringiro

Nicolas Andre↵
email: nicolas.andre↵@femto-st.fr

Calibration of a multi-mobile coil magnetic manipulation system utilizing a control-oriented magnetic model

Keywords: calibration, multi-mobile coil, magnetic manipulation, microrobotics, numerical algorithm

published or not. The documents may come

Calibration of a multi-mobile coil magnetic manipulation system utilizing a control-oriented magnetic model

Nahashon Osinde, Maxime Etievant, Jean Bosco Byiringiro, Nicolas Andreff

Introduction

In the last decade, there has been noteworthy development in the field of microrobotics, largely attributed to the advancements of micro-fabrication techniques and microscopy [START_REF]Microbiorobotics: biologically inspired microscale robotic systems[END_REF]. Researchers are exploring various ways and techniques in which systems can be developed to perform microscale tasks. Microrobots have several applications in engineering as well as in the medical field. These include micro-manipulation, microassembly, performing minimally invasive surgical operations and performing targeted diagnosis and drug delivery [START_REF] Sikorski | The ARMM System: An Optimized Mobile Electromagnetic Coil for Non-Linear Actuation of Flexible Surgical Instruments[END_REF][START_REF] Abbott | Magnetic Methods in Robotics, Annual Review of Control[END_REF][START_REF] Nelson | Microrobots for minimally invasive medicine[END_REF]. Applications such as minimally invasive surgery have several benefits which include; reduction in recovery time, medical complication and infection risks [START_REF] Mack | Minimally Invasive and Robotic Surgery[END_REF].

With the size of microrobots, there is di culty in controlling them due to viscous forces dominating inertia forces at the micro-scale. Therefore, building of robots in the microscale has technical limitations to engineers since scaling down of concepts and designs used in the macro-scale is di cult and at times impractical. Engineers therefore have to abandon intuition gained from designing in the macro-world and instead rely on the physics governing microrobot design and analysis. Descending in scale, it becomes more and more interesting to eliminate energy storage elements and introduce remote actuation rather than design a sophisticated microrobot. An e cient way to control microrobots wirelessly is by use of magnetic fields. In addition, magnetic fields are easy to produce, have limited interaction with the human body and are able to transfer large amounts of power [START_REF] Sikorski | Vision-Based 3-D Control of Magnetically Actuated Catheter Using BigMag-An Array of Mobile Electromagnetic Coils[END_REF].

Electromagnetic coils are a suitable choice in generating magnetic fields since they are inactive when powered down [START_REF] Yang | DeltaMag: An Electromagnetic Manipulation System with Parallel Mobile Coils[END_REF]. By modulating the magnetic field, it is possible to generate a force and torque which can be used in the control of a microrobot. In particular, higher dexterity and e ciency can be achieved by modulating the magnetic field through both current variation and motion of the coils as initially shown in [START_REF] Veron | Dealing with redundancy of a multiple mobile coil magnetic manipulator: A 3rpr magnetic parallel kinematics manipulator[END_REF] and later exploited in [START_REF] Sikorski | Introducing BigMag -A novel system for 3D magnetic actuation of flexible surgical manipulators[END_REF][START_REF] Chen | Study of robotized electromagnetic actuation system for magnetic microrobots devoted to minimally invasive ophthalmic surgery[END_REF][START_REF] Chen | Mathematical approach for the design configuration of magnetic system with multiple electromagnets[END_REF]. The use of mobile coils present several advantages over the use of stationary bigger coils when a large footprint of the patient workspace is needed. Large coils which are heavy may result in complexities in the structural design thereby making the manufacturing process more expensive [START_REF] Chautems | The tethered magnet: Force and 5-DOF pose control for cardiac ablation[END_REF]. In addition, larger coils would have higher inductance which would limit the control frequency, increase the power consumption and the heat dissipated. Therefore, in a bid to have a large workspace with minimal energy expenditure, it is essential to use mobile coils and move them close to the area of manipulation.

Consider a magnetic manipulation system having several coils made by mobile serial robots as seen in Figure 1. For highly accurate positioning and control of a microrobot through magnetic manipulation, a control-oriented model that is fast and accurate to compute the magnetic field is required [START_REF] Chowdhury | Controlling multiple microrobots: recent progress and future challenges[END_REF]. Accurate modeling of the magnetic field can be achieved using Finite Element Models but this is too slow for control. For control, one usually relies either on the simplified dipole model or on prerecorded magnetic field. It has been shown in [START_REF] Petruska | Optimal permanent-magnet geometries for dipole field approximation[END_REF] that the dipole model is accurate for distances far away from the magnet. These approaches are not perfectly suited to multi-mobile electromagnetic manipulation. To perform multi-mobile electromagnetic manipulation closed-loop control, a control-oriented model is needed, i.e. a model that provides a trade-o↵ between the computation time and accuracy. One such model based on complete elliptic integrals has been shown in [START_REF] Etievant | An Improved Control-Oriented Modeling of the Magnetic Field[END_REF] to be suitable for closed-loop control of magnetically actuated robots. A cost e↵ective way to improve this model further is through modelbased calibration [START_REF]Robot calibration, 1st Edition[END_REF][START_REF] Johansson | System Modeling and Identification, Information and system sciences series[END_REF]. Calibration is a process through which parameters to a given system model are identified and adjusted such that the output of the model fits well to the actual output of the system [START_REF] Elatta | An Overview of Robot Calibration[END_REF]. Through calibration of an electromagnetic system, the positioning accuracy of a microrobot can thus be enhanced by adjusting model parameters rather than altering the mechanical structure of the electromagnetic system. In [START_REF] Petruska | Model-Based Calibration for Magnetic Manipulation[END_REF] a generic magnetic model based on a multi-source spherical multipole expansion was calibrated using nonlinear least squares. Although numerically e cient, the model used in [START_REF] Petruska | Model-Based Calibration for Magnetic Manipulation[END_REF] was less grounded into physics as the model in [START_REF] Etievant | An Improved Control-Oriented Modeling of the Magnetic Field[END_REF] and was also not assessed in close vicinity of the coils.

In this paper a model of the full multi-mobile magnetic manipulation system (including kinematics of the manipulator) is presented. It is shown that the calibration of the system needs to be split into several subsystems which eventually results in calibrating the coil and a reference rigid transformation only. The coil was modelled using the control-oriented model which provides a trade-o↵ between its accuracy and computation speed [START_REF] Etievant | An Improved Control-Oriented Modeling of the Magnetic Field[END_REF]. Sensitivity analysis on the model to investigate the significance of the model parameters is also presented.

This paper is composed of six sections. We first present the background on calibration in Section 2. Modelling of the magnetic manipulation system is highlighted in Section 3. Section 4 presents calibration of the mobile magnetic manipulation system and Section 5 is devoted to the calibration results. Conclusions and perspectives are drawn in Section 6.

Calibration background

Calibration has several phases which include; modelling, measurement, parameter identification, error compensation and finally validation [START_REF] Ljung | System identification: theory for the user, 2nd Edition, Prentice Hall information and system sciences series[END_REF]. The steps of calibration are described briefly in the subsequent sub-topics with specific insight on magnetic manipulation.

Modelling

The first step of calibration in this context is the formulation of a mathematical model that describes the magnetic field. Several models are in use for the computation of the magnetic field [START_REF] Chen | Mathematical approach for the design configuration of magnetic system with multiple electromagnets[END_REF]. Nonetheless, for real time control a model that is both accurate and fast to compute is needed. It has been shown in [START_REF] Etievant | An Improved Control-Oriented Modeling of the Magnetic Field[END_REF] that the magnetic model based on complete elliptic integrals is accurate in the whole workspace area and takes less time to compute as compared to existing magnetic models. Moreover, when the coils are mobile, a kinematic model of the moving structure is needed.

Measurement

This step involves collection of magnetic field measurement data that will be used in the calibration process. This data is obtained by measuring the magnetic field generated by the coil in the workspace using a Hall e↵ect sensor.

Parameter identification

Parameter identification can be defined as the numerical procedure used to determine parameter values to the magnetic model. This step has been addressed by several researchers and various parameter identification algorithms have been used. In this paper, we rely on the output error identification method [START_REF] Ljung | System identification: theory for the user, 2nd Edition, Prentice Hall information and system sciences series[END_REF][START_REF] Söderström | System Identification[END_REF] as shown in Figure 2.

Let ⇠ define the set of parameters in the model, the aim is to minimize the objective function generally defined from the root mean square error (RMSE). The RMSE highlights the discrepancy between the magnetic field computed from the model and the one obtained from the measurement. To compute the deviation, ✏ i , at a particular measurement configuration, i, in the workspace, equation (1) is used:

✏ i (⇠) = Bi (I, w T c , ⇠) B i (I, w T c ) ( 1
)
where Bi is the estimated magnetic field from the model using the control inputs (I the current in the coil, w T c the pose of the coil in space) and the model parameters, ⇠. B i is the measured output. RMSE from several m points in the workspace is therefore computed using equation [START_REF] Sikorski | The ARMM System: An Optimized Mobile Electromagnetic Coil for Non-Linear Actuation of Flexible Surgical Instruments[END_REF]:

Multi

RMS E = 1 2 v t 1 m m X i=1 k✏ i k 2 (2)
The resulting cost function that needs to be minimized using a numerical algorithm so as to obtain the optimal parameter values ⇠ ⇤ is:

⇠ ⇤ = argmin {⇠} 1 2 v t 1 m m X i=1 k✏ i k 2 (3)
Algorithms used in parameter identification include Sequential Quadratic Programming (SQP), Nelder-Mead, Levenberg-Marquardt, genetic algorithm among many others [START_REF] Nocedal | Numerical Optimization[END_REF].

Error compensation

Error compensation (also referred to as the correction step) involves adjusting the parameter values in the control-oriented model with the values obtained after parameter identification. This is demonstrated by the arrow crossing the 'Model' box in Figure 2.

Validation

This is the final step in the calibration process. It involves comparison of the magnetic field obtained from the modified model and the measured one. Measurement data used in this step must be statistically independent from the data used in parameter identification. The error computed from equation ( 2) should be as minimum as possible, so as to ascertain that the model is valid. If the model is not validated, parameter identification is done again, or new measurement data are taken or a new model of the electromagnetic system is formulated again. A generic case study will be first developed. A schematic of a magnetic manipulation system is as shown in Figure 3. It consists of several serial robots, each with a magnetic coil attached to the end-e↵ector. Further developing this generic case study, each coil generates a magnetic field B j . The total field at a position configuration P in the world reference frame for n coils, assuming the superposition theorem holds is w B( w P) = P n j=1 w B j ( w P) [START_REF] Abbott | Magnetic Methods in Robotics, Annual Review of Control[END_REF]. To obtain B j , modelling of one circular coil loop, the coil and the robotic manipulator needs to be done.

Modelling of the magnetic manipulation system

Generic case

Modelling of one circular coil loop

Consider a circular current loop of one coil in the spherical coordinate system with loop radius a centered at O as shown in Figure 4. The magnetic field in a point P of spherical coordinates (r, ✓, ) computed using complete elliptic integrals is [START_REF] Etievant | An Improved Control-Oriented Modeling of the Magnetic Field[END_REF]:

8 > > > > > > > > > > < > > > > > > > > > > : B r (r, ✓) = µ 0 I ⇡ a 2 p a 2 +r 2 +2arsin(✓) E(k)cos(✓) a 2 +r 2 2arsin(✓) B ✓ (r, ✓) = µ 0 I ⇡ 1 p a 2 +r 2 +2arsin(✓) f (r, ✓) B = 0 f (r, ✓) = [ E(k)(r 2 +a 2 cos(2✓)) 2sin(✓)(a 2 +r 2 2arsin(✓)) K(k) 2sin(✓) ] (4)
where µ 0 is the vacuum permittivity, K(k) and E(k) are the complete elliptic integral functions of first and second kind:

K(k) = Z ⇡ 2 0 1 p 1 k 2 sin 2 ✓ d✓ E(k) = Z ⇡ 2 0 p 1 k 2 sin 2 ✓d✓ ( 5
)
and k is the modulus of the elliptic integral functions:

k 2 = 4arsin✓ a 2 + r 2 + 2arsin✓ (6)
However, the expressions in equation ( 4) are not defined for k = 0. When k << 1, the continuity is done by using [START_REF] Schill | General relation for the vector magnetic field of a circular current loop: a closer look[END_REF]:

K(k) = ⇡ 2 + k 2 ⇡ 8 + k 4 9⇡ 128 + 0(k 4 ) E(k) = ⇡ 2 k 2 ⇡ 8 k 4 3⇡ 128 + 0(k 4 ) (7)
With the continuity expression in equation ( 7), equation ( 4) is simplified and the magnetic field for k << 1 becomes:

8 > > > > > > < > > > > > > : B r (r, ✓) = µ 0 Ia 2 4 cos(✓) (2a 2 +sin(✓)ar+2r 2 ) (a 2 +2sin(✓)ar+r 2 ) 5 2 B ✓ (r, ✓) = µ 0 Ia 2 4 sin(✓) (2a 2 +sin(✓)ar r 2 ) (a 2 +2sin(✓)ar+r 2 ) 5 2 B = 0 (8)
In addition, note for consistency sake that when r >> a, the expression in equation ( 8) is simplified into the dipole approximation [START_REF] Grant | Electromagnetism, 2nd Edition, The Manchester physics series[END_REF].

Modelling of one coil

The magnetic coil is composed of many single loops. Figure 5 shows the actual coil and its schematic used in this case study. Although the actual wiring is not perfectly known or made, we assume a regular model as shown in Figure 6, where d is the diameter of the copper wire used in the winding. Copper wire is wound cylindrically as seen in Figure 5a. Loops of wire wound in the radial direction are referred to as radial loops and loops of wire along the axis are referred to as axial loops. They are shown in the schematic of the coil in Figure 5b. For modelling and later calibration to be done on the coil, information from the model given in equation ( 4) needs to be extracted [START_REF] Renaud | Simplifying the kinematic calibration of parallel mechanisms using vision-based metrology[END_REF]. These comprises of: constants, measured variables and parameters to be estimated. From the model, the constant is vacuum permittivity, µ, measured variables are the po- sition of the measurement point (obtained from r and ✓) and the parameters to be estimated are:

?d = 2R w R c N r ⇥ d N a ⇥ d Z Y X
• External loop radius (R c ).

• Wire radius (R w ).

• Number of loops radially (N r ).

• Number of loops axially (N a ).

These parameters need to be taken into account in equation ( 4) so as to compute the total magnetic field in the coil. The parameter values to the assumed regular model of the coil might not be correct, hence the need for calibration to identify them.

Expressing the radius, a and height, dz of each loop as a function of its index in the radial wiring, n r and axial wiring,

n a : a(n r ) = (R c R w ) (2R w (n r 1)) (9)
dz(n a ) = (2R w (n a 1)) + R w [START_REF] Chen | Study of robotized electromagnetic actuation system for magnetic microrobots devoted to minimally invasive ophthalmic surgery[END_REF] There is also need to compute r and ✓ for each loop in the coil:

r(P, n a ) = (X 2 + Y 2 + (Z + dz(n a )) 2 ) 1 2 (11) ✓(P, n a ) = arccos( Z + dz(n a ) r(P, n a ) ) (12)
where P = (X, Y, Z) T are the coordinates of the point P in the reference frame F c attached to the coil. The magnetic field sum for each loop axially and radially in the coil reference frame is then considered to compute the total magnetic field:

8 > > > > > > > > < > > > > > > > > :

Modelling of the full multi-mobile magnetic manipulation system

The multi-mobile electromagnetic manipulation system is first subdivided into several arms. Figure 7 shows the schematic representation of one arm in two configurations. F w is the world coordinate system attached to the static Hall e↵ect sensor used for calibration, F b j the j'th robot base frame, F e j the j'th robot end-e↵ector frame and F c j is the j'th robot coil frame. Expressing the magnetic field in the world frame at the origin of F w (measurement point):

w B j ( w Õ) = w T c j c j B j (⇠ m, j , w T 1 c j w Õ) (14)
where B = B 0

! , w Õ = w O 1
! are the projective coordinates of vector B and point O , ⇠ m, j is the vector of coil model parameters and w T c j = w R c j w t c j 0 1 ! which is the pose of the coil frame relative to the world frame given by the kinematic chain:

w T c j = w T b j FKM(q j , ⇠ g, j ) e j T c j ( 15
)
where FKM is the forward kinematic model of the robot, q j , the variable robot joint values and ⇠ g, j , the constant kinematic parameters to the robot. w T b j and e j T c j are also constant transformations.

The total magnetic field at the origin of the Hall e↵ect sensor in the world frame from all the arms of the multi-mobile electromagnetic system is therefore given by:

w B( w T b j , q j , ⇠ g, j , e j T c j , ⇠ m, j , w Õ) = n X j=1 w T c j (•) c j B j (⇠ m, j , w T c j (•) 1 w Õ) ( 16
)
where w T c j (•) means that equation ( 15) is used to compute w T c j using the variable q j and constant parameters w T b j , e j T c j and ⇠ g, j .

Calibration of the mobile magnetic manipulation system

Using equation ( 16) directly in the calibration cost function would yield a high dimensional problem. Careful examination of the equation that the numerical problem can be simplified. In building model, we assumed as most researchers [3] that the inter-coil interferences are neglectable and that the superposition theorem holds. Additionally, since the coils can be individually turned on and o↵, a reasonable calibration protocol is to turn them on one at a time. As a consequence, only one term in equation ( 16) will be non-zero at a time and each arm can thus be simplified in n subcalibration problems for each m configuration.

8 j = 1..n argmin { w T b j , e j T c j ,⇠ g, j ,⇠ m, j } 1 2 v t 1 m m X i=1 k w B j (i) w T c j (•) (i) c j B j (i) (•)k 2 (17)

Robot forward kinematic calibration

To use the FKM, the robot needs to be calibrated. Calibration of the serial robot necessitates moving it in its whole workspace and taking measurements. Using directly equation ( 16), a full calibration of both kinematic and magnetic parts of the model is not possible. Indeed, the coils mounted on the end-e↵ector generate a magnetic field in a limited workspace (near the coil). Measurements of the magnetic field in the serial robot workspace which is significantly larger than the coil workspace would be zero at most locations. To get relevant data on the magnetic field, we are therefore restricted to a small workspace i.e. near the coil. On that account, the calibration of the kinematic model of the robot needs to be done separately.

Robot-world and hand-coil calibration

With the assumption that the robot is already calibrated, what remains of importance is the calibration of the coil as well as the robot base to world ( w T b j ) and hand-coil ( e j T c j ) transformations. Determining these two transformations is a problem similar to the simultaneous robot-world and hand-eye calibration problem in [START_REF] Horaud | Hand-eye calibration[END_REF]. However, since we cannot compute or measure the full 4 ⇥ 4 transformation w T c j but only the 3 ⇥ 1 magnetic field value at the sensor point B( w T c j) , we need to adapt the calibration method to this partial information through further simplification.

Simplified calibration problem

Taking into account all the considerations above, the remaining unknown parameters are the transformation w T c j of the j'th robot coil with respect to the sensor and ⇠ m, j the magnetic parameters to the j'th coil model. Keeping the idea of hand-eye calibration [START_REF] Andre↵ | Robot hand-eye calibration using structure-from-motion[END_REF], a 'relative' model using relative displacement of the coil with respect to a reference configuration instead of an 'absolute' model is written. Consider configuration (0) and (i) after motion of the robot arm in Figure 7, the relative model is given by expressing w T c j (i) (•) as:

w T c j (i) (•) = w T c j (0) e j T c j 1 b j
T e j (q j (0) , ⇠ g, j ) 1 b j T e j (q j (i) , ⇠ g, j ) e j T c j (18) Thus,

w T c j (i) (•) = w T c j (0) e j T c j 1 (0) T (i) e j T c j (19)
where (0) T (i) is the displacement of the end-e↵ector from F e j (0)

to F e j (i) . For the experimental protocol, we propose to use only pure translations such that R = I, therefore, (0) T (i) = I (0) t (i) 0 1

! . Indeed, it makes it simpler for the user and reduces again the complexity of the calibration. Under this protocol, equation ( 19) develops as:

w T c j di f f (•) = w T c j (0) I e j R c j T (0) t (i) 0 1 ! ( 20
)
and thus:

w T c j di f f (•) 1 = I e j R c j T t 0 1 ! w T c j (0) 1 = w R c j T (0) e j R c j T (0) t (i) w R c j T (1) w t c j (0) 0 1 ! (21)
Since we have the computer aided design (CAD) values of the mounting brackets of the coil to the end-e↵ector, we can reasonably estimate the hand-coil transformation up to an accuracy which is in the order of magnitude of the FKM accuracy. Depending on the machine used in the fabrication of the mounting brackets, the hand-coil transformation accuracy can even be better.

The error on e j R c j T (0) t (i) is therefore of the same order of magnitude as the error on t. It can be noticed that e j t c j vanishes and only e j R c j is needed. We can then use the CAD values for e j R c j T . Moreover, in most set-ups, e j R c j will be a ⇡ 2 rotation or the identity. Consequently, removing the projective notation the calibration becomes:

argmin { w T c j (0) ,⇠ m } 1 2 v u u u u u u u u t 1 m m X i=1 k w B(i) w R c j (0)c j B j (i) (⇠ m, j , ( e j R c j T (0) t (i) w R c j (0) T w t c j (0) ))k 2 ( 22
)
where w B(i) is the measured magnetic field at configuration (i) of the robot using the Hall e↵ect sensor and c j B j (i) is the estimated magnetic field at configuration (i) using the model. The values to w T c j (0) and ⇠ m are estimated through calibration.

Sensitivity analysis of the coil parameters

The coil model has many parameters and not all of them may have a major influence on the output variability. In a bid to understand the significance of each parameter, a sensitivity analysis was conducted. This is a fundamental step before calibration [START_REF] Ray | Modelling nitrogen and carbon cycles in Hooghly estuary along with adjacent mangrove ecosystem[END_REF]. The elemental approach to performing sensitivity analysis is utilizing partial di↵erentiation and the simplest one, varying one factor at a time (OAT). Depending on the nature of the model, partial di↵erentiation may become so complex rendering implementation impractical. OAT is referred to as 'local sensitivity' since the analysis is only done for the selected point estimates and not the entire distribution space. Another approach to sensitivity analysis is through factorial designs. Using factorial designs for sensitivity analysis involves choosing a parameter distribution space and implementing the model for all combinations in the distribution space. However, as the number of parameters and size of the distribution space increases, so does the number of runs. Nevertheless, the same information can be obtained through fractional factorial design approaches such as Box-Behnken and Taguchi [START_REF] Robles | Design of experiments for sensitivity analysis of a hydrogen supply chain design model[END_REF]. For this research, factorial design was used, it has the following advantages:

• Parameter interaction e↵ects are able to be analysed and assessed. This leads to good understanding of how the output response changes with the input parameters.

• Output response analysis is able to be done at di↵erent parameter levels yielding valid conclusions over a wide range of experimental conditions.

A full factorial design with seven levels of unitary step for the parameter distribution space as shown in It was done using simulation data based on the complete elliptic integral model given in equation [START_REF] Chowdhury | Controlling multiple microrobots: recent progress and future challenges[END_REF]. To draw conclusions from the sensitivity analysis, main e↵ects and interaction e↵ects plots were studied. The main e↵ects plot given in Figure 8a describes how the RMSE, defined in equation ( 2) changes as the value of one parameter is varied. The x-axis on the plot highlights the distribution space through which each parameter is swept across.

How parallel the graph of a parameter is to the x-axis, indicates the degree of variability in the model's output due to the considered distribution space of that parameter. The steeper the slope, the greater the degree of variability due to that input parameter. As seen from the plot, the external radius has the greatest main e↵ect and the wire radius has the least main e↵ect. However, to make a valid conclusion on the most influential parameter, the interaction e↵ects between the parameters need to be analysed. The interaction e↵ects plot is shown in Figure 8b.

In examining the interaction e↵ects plot, lines for two parameters that appear to intersect indicate presence of interaction while lines that appear to be translated and not overlap indicate no interaction between the parameters. It can be observed that, there is no interaction between the external radius and number of loops radially. However, there is a strong interaction between the external radius and wire radius and the number of loops axially and radially. The strong interaction between the external radius and wire radius can be attributed to the fact that the equation used in obtaining the external radius is a function of the wire radius. Hence, only one of them can be used in practice. The model can be simplified by excluding the wire radius because its value can be obtained from measurement data of the manufacturer.

Numerical solver selection

For parameter identification, several algorithms and optimization methods have been put forth. These algorithms could be categorized into; gradient methods, gradient-free methods and genetic based methods [START_REF] Lo-Amni | Automated Calibration for Numerical Models of Riverflow[END_REF]. One algorithm from each of these categories was selected for this case study. Sequential quadratic programming (SQP), Nelder-Mead algorithm and genetic algorithm, all implemented in MATLAB 2019a were considered for calibration. SQP is an e↵ective method used in the optimization of small and large problems having pronounced non-linearities [START_REF] Nocedal | Numerical Optimization[END_REF]. Nelder-Mead algorithm has been used in parameter estimation where the values of the function are uncertain and subject to noise [START_REF] Larson | Derivative-free optimization methods[END_REF]. Genetic algorithm has been used to solve optimization problems in which the objective function is discontinuous, stochastic, non-di↵erentiable or nonlinear [START_REF] Tao | Configurable Intelligent Optimization Algorithm[END_REF]. To analyse the performance of these algorithms, the metric used was the root mean square error (RMSE).

To select which parameter identification method to use, calibration was first performed using the three di↵erent algorithms using simulated data; sequential quadratic programming (SQP) implemented as fmincon in MATLAB, Nelder-Mead algorithm implemented as fminsearch in MATLAB and genetic algorithm. Calibration was done using simulation data on a 5⇥5⇥5 regular grid.

Simulation data was obtained numerically from the model based on complete elliptic integrals in order to assess the convergence properties of each algorithm. The starting values in implementing the algorithms were the same for all the solvers. Table 2 shows the results of the RMSE percent reduction after calibration using each of the numerical algorithms.

SQP Genetic algorithm

Nelder-Mead algorithm RMSE % reduction 88.4830% 94.1784% 96.4259% It is observed from Table 2 that Nelder-Mead algorithm has the best performance when comparison is done using the percentage reduction of the RMSE. This is followed by genetic algorithm and finally SQP. The performance of SQP could be attributed to using the numerical Jacobian. A lower RMSE indicates that the magnetic field computed by the model closely matches the measured magnetic field. It is desirable that the RMSE is as small as possible. In every calibration done using simulation data, the Nelder-Mead algorithm outperformed the other two. Nelder-Mead algorithm is an unconstrained direct search algorithm. Unlike more traditional optimization methods that use information about the gradient or higher derivatives to search for an optimal point, a direct search algorithm searches a set of points around the current point, looking for one where the value of the objective function is lower than the value at the current point. Direct search methods are used to solve problems for which the objective function is not di↵erentiable, or is discontinuous. The Nelder-Mead method has a recommendable performance in the vicinity of the local or global minimum since the intensive computation of the first and second derivative can be avoided [START_REF] Lo-Amni | Automated Calibration for Numerical Models of Riverflow[END_REF].

Calibration results

Simulation of the coil calibration

For simulation purposes the transformation matrix w T c j

in equation ( 22) was dropped and calibration done only on the parameters to the coil model, ⇠ m . Data used for the simulation was obtained numerically from the model based on complete elliptic integrals proposed in [START_REF] Etievant | An Improved Control-Oriented Modeling of the Magnetic Field[END_REF] which was shown to be close and accurate enough in computing the magnetic field. Validation was done using simulated data on a 20 ⇥ 20 ⇥ 20 regular grid with similar dimensions as the one used for calibration. The ground-truth values, the parameter values before (values used to initialize the numerical algorithm) and after calibration are given in Table 3. error is the deviation divided by the measured magnetic field at a particular point in the workspace. A higher relative error value indicates a huge deviation at that particular point. From Figure 9a, it is seen that there is a huge magnitude of relative errors on the sides of the coil as opposed to the rest of the workspace.

Ground

Regardless of the huge errors on the side of the coil, it is not a point of concern since the microrobot during magnetic manipulation will not be placed at that location Plots of the angular error in the XZ and the XY plane are as shown in Figure 10a and Figure 10b respectively. The angular deviation illustrates the error between the orientation of the computed field and the measured one. An angular deviation induces an error on both the computed force and torque exerted on the robot during manipulation. It is preferable to have the angular deviation as low as possible. From the colourmap plot in Figure 10a, a comparison can be drawn to the colourmap plot of relative errors discussed in an earlier paragraph. Majority of the huge magnitude of errors are on the sides of the coil as opposed to the interested workspace. As has been discussed, the huge errors on the side are not a point of concern during magnetic manipulation since the robot will not be placed at that particular location.

Experiment

For the physical experiment, magnetic field measurement data was collected in a 5 ⇥ 5 ⇥ 5 as well as in a 8 ⇥ 8 ⇥ 8 equally spaced grid. To generate the measurement grid, the coil was translated in its X, Y, Z direction with the help of a serial robot

Conclusions and perspectives

A calibration protocol used to calibrate a magnetic manipulation system with multiple mobile coils has been presented. The magnetic manipulation system was subdivided and each coil calibrated separately assuming the superposition theorem holds and that the robotic arms are perfectly calibrated. A controloriented magnetic model based on complete elliptic integrals with its continuity extension was used. The protocol involved identifying model parameter values to the electromagnetic system by use of the Nelder-Mead algorithm. Calibration was performed on both numerically generated data, as well as data measured from a physical system. It was observed that the root mean square error was reduced after calibration indicating an increase on the model accuracy. A minimum as possible value of the root mean square error is desired. In this research, it was assumed that the kinematic parameters to the serial robot were perfectly known. In a case where they are not known, the serial robot needs to be decoupled from the coil and its calibration done separately. It was also assumed that the superposition theorem holds and that there are no couplings between the electromagnets. This assumption might not always be true especially when the coils are too close to each other and could open up further research paths related to electromagnetics. In future works, advanced control strategies to magnetically navigate surgical instruments in relevant medical procedures will be developed. The capabilities of our magnetic manipulation system will also be improved by introducing ferromagnetic-cores in the coils so as to achieve stronger magnetic fields. The introduction of such cores would necessitate modification of our control-oriented model.

Figure 1 :

 1 Figure 1: Magnetic manipulation of a microrobot using two coils. Each coil is being displaced by a UR3 (Universal Robot) serial robot.

Figure 2 :

 2 Figure 2: Calibration process through output error minimisation.

Figure 3 :

 3 Figure 3: Schematic of a magnetic manipulation system.

Figure 4 :

 4 Figure 4: Circular current coil loop.

aFigure 5 :

 5 Figure 5: (5a) Side view of the actual coil. (5b) Schematic representation of the various parameters that need to be identified during the calibration phase.

Figure 6 :

 6 Figure 6: Side view of a regular model of the coil.

  r(P, n a ), ✓(P, n a ), a(n r )) c B ✓ (P) = Na P na=1 Nr P nr=1 B ✓ (r(P, n a ), ✓(P, n a ), a(n r )) c B (P) = 0 (13)

Figure 7 :

 7 Figure 7: Schematic of one arm of the magnetic manipulation system in a reference (0) configuration and in a generic (i) one.

4. 1 .

 1 Simplifying the calibration problem 4.1.1. Separate calibration of each arm

Figure 8 :

 8 Figure 8: (8a) Main e↵ects plot. How parallel the graph of a parameter is to the x-axis, indicates the degree of variability in the model's output due to the considered distribution space of that parameter. (8b) Interaction e↵ects plot. Lines for two parameters that appear to intersect indicate presence of interaction while lines that appear to be translated and not overlap indicate no interaction between the parameters.

Figure 10 :

 10 Figure 10: Angular errors after calibration using simulated data.

Figure 12 :

 12 Figure 12: Relative errors after calibration using experimental data.

Figure 13 :

 13 Figure 13: Angular errors after calibration using experimental data.

Table 1 :

 1 Table 1 was done using Minitab 17 Statistical software. Parameter distribution space used in factorial design.

	Parameter	Parameter distribution space
	External loop radius (R c )	48E 03m to 54E 03m
	Wire radius (R w )	46E 05m to 52E 05m
	Number of loops axially (N a )	45 to 51
	Number of loops radially (N r )	22 to 28

Table 2 :

 2 Algorithm comparison using simulated data.

Table 3: Ground-truth values, parameter values before and after calibration using simulation data.The following statistical metrics were then studied; maximum error, RMSE and relative error. For a better cross-validation, an additional metric was used, the angular deviation between the measured and modelled fields. The angular deviation is independent from the metrics used for the identification, both by the data it uses and by it's mere definition. Table4shows the errors obtained after calibration while Table5shows the percent reduction. As expected, the errors after calibration are less than before calibration indicating an improvement in the parameter estimation of the magnetic model.

Acknowledgements

This work is part of Multiflag project (ANR-16-CE33-0019) funded by the French National Research Agency (ANR). It was also supported by EIPHI Graduate School (Contract No. ANR-17-EURE-0002), by Région Bourgogne Franche-Comté through the COErCIVe project and by the Grand Prix Scientifique 2018, Fondation Charles De↵orey, Institut de France.

Absolute error

Before calibration After calibration Maximum error(T) 0.0049 0.0033 Root mean square error(T) 1.1000E 03 4.2622E 04 Maximum error is the maximum deviation between the computed magnetic field and the measured magnetic field. A typical location of where the maximum error might be pronounced is on the sides of the coil as opposed to the interested workspace whose dimension is 100mm by 100mm by 50mm. This can be observed in Figure 9.

Table 6 shows the errors obtained in the interested workspace after calibration while Table 7 shows the percent reduction. A small RMSE indicates a tight fit of the model to the measurement data. A small as possible RMSE is desired. It is observed that after calibration the RMSE was reduced by 65%. To get a more precise look of where the deviations are, a colourmap of the relative errors and angular deviation in the workspace was plotted. Plots of the relative error in the XZ and the XY plane are as shown in Figure 9a and Figure 9b respectively. Relative (Universal robot 3) arm while the sensor was held in a fixed position on the worktable as shown in Figure 11. The former grid was used for calibration while the latter for validation. The data was collected in a workspace of length, 100mm, width 100mm and height, 50mm. This workspace is close to the coil and is the area of interest in magnetic manipulation. Measurement data was collected using a Hall e↵ect sensor, THM1176 from Metrolab. It is a 3-axis Hall magnetometer that has a range from 100µT to 14T . The sensor also simultaneously measured the magnetic field in all the three axes (X, Y, Z). Current was provided to the coil through a power supply and to ascertain that the right current magnitude was being delivered to the coil, a multimeter was connected in the circuit. EZMag3D software from Metrolab was then used in data acquisition of the magnetic field from the sensor.

Absolute error

Parameter identification using equation ( 22) was done in MATLAB using the Nelder-Mead algorithm. Parameters of the magnetic coil, ⇠ m and the transformation matrix w T c j (0) obtained from physical measurements used to initialize the algorithm and the values after calibration are given in Table 8. Similarly to the simulation results, Table 9 shows the errors obtained after calibration while Table 10 shows the percent reduction. On performing validation, it was observed that the RMSE between the model and measurements was reduced by 37% after calibration.

Absolute error

Before calibration After calibration Maximum error(T) 0.0092 0.0074 Root mean square error(T) 0.0043 0.0027 Plot of the relative error and the angular deviation in the interested workspace are as shown in Figure 12 and Figure 13 respectively. The magnitude of errors obtained using experimental data was observed to be much greater than those obtained using simulation data. This may be attributed to the positioning accuracy of the serial robot and numerical errors that may have arisen during processing of the collected data.