3.1. Numerical analysis of symmetric and asymmetric configurations

To understand how positioning the electromagnets affects the induced magnetic field, two sets of numerical analysis are performed for the symmetric and asymmetric configurations. The existence of the ferromagnetic cores and the large air gap between the two electromagnets complicate the analysis of the generated magnetic field. Therefore, calculations were carried out with the aid of an appropriate PDE solver program.

The numerical model analysis was built based on an optimized electromagnet design presented in a previous study of us for the purpose of the proposed EMA system [Reference Alasli, Çetin, Akçura, Kahveci, Can and Tamer35]. The electromagnet was enforced with a ferromagnetic core in order to first increase the intensity of the magnetic field since the proposed implementation is at micro scale and second to enlarge the untethered workspace of the microrobot. Moreover, the core of the electromagnets was specifically designed to minimize the inhomogeneity of the magnetic field at the symmetric configuration and the gradient deviation around the central axis at the asymmetric configuration. The current of electromagnets is kept constant (as 2 A) for all analysis since the effect of the positioning is investigated in this study.

3.1.1. Symmetric configuration

To analyze the induced electromagnetic field at the symmetric configuration, a parametric simulation is carried out at three different distances $D=D_1=D_2$ of 50, 60, and 70 mm. The obtained flux intensities $\vec{B}$ on the coaxis of the pair as function of the distance $D$ are plotted in Fig. 4(a). As can be seen, the resultant magnetic field is a superposition of the two magnetic fields arising from each electromagnet. It can also be observed that with an increment in the distance between the electromagnets, a decrement in the gradient near the center point is accompanied by a decrement in the intensity of the magnetic field. It can also be noted that it is not possible to obtain a perfect zero gradient region with the proposed configuration as any other coil pair setup. Therefore, a low gradient region is assumed to be a homogeneous field region if the induced force due to the gradient at any point of this region is lower than the friction force acting on the microrobot. From the aspect of homogeneity, if we consider the maximum accepted deviation in the gradient on the central line to be 5%, the width of the corresponding homogeneous line segments on the axis is measured as 16.4, 17.6, and 20.1 mm at the symmetric configuration for $D_1=D_2=D$ equals to 50, 60, and 70 mm, respectively. The results also show that this quasi-homogeneous region at the center can be widened by increasing the distance $D$ . Furthermore, it can be realized that the symmetric configuration can be used as an alternative to Helmholtz Coil due to its advantage of increasing the intensity of magnetic field and having an adjustable workspace.

Figure 4.

Magnetic flux density on central axis in (a) symmetric configuration and (b) asymmetric configuration at different distances.

3.1.2. Asymmetric configuration

A similar parametric analysis is performed for the asymmetric configuration. The distances from the center point to each electromagnet are defined as $D_1=D$ and $D_2=D+\epsilon$ . The distance $D$ is set to 60 mm, and the asymmetry parameter ( $\epsilon$ ) is set to 0, 10, 20, 30, and 40 mm ( $D_2=$ 60, 70, 80, 90, 100 mm, respectively). The obtained intensity of the magnetic field $\vec{B}$ on the central axis is given as a function of the distances, as shown in Fig. 4(b). As it can be realized, the form of the field varies with the asymmetric positioning of the electromagnets ( $\epsilon$ ) when it is compared to the symmetric case. Moreover, the point at which the minimum magnetic field intensity occurred is shifting toward the farther electromagnet. Hence, it can be stated that the asymmetric placement of the electromagnets can be used to create a magnetic field gradient similar to Maxwell coil in any predefined region.

Furthermore, the two previous simulations’ results confirm that by using this re-configurable coaxial electromagnets set, the magnitude and gradient of the magnetic field can be modified purposefully. Namely, the proposed setup has the capability to create a combined effect of Maxwell and Helmholtz coils.

Figure 5.

Magnetic field gradient at center point $G_u(\epsilon )$ with respect to asymmetry parameter $\epsilon$ (numerical results).

3.2. Relationship between the asymmetry parameter and the magnitude of magnetic force

The slope of the axial magnetic field gradient ( $\frac{\partial B}{\partial u}$ ) at the center point $O$ in Fig. 4(b) increases while the pair of electromagnets transforms from symmetrical to asymmetrical configuration (with $\epsilon$ ). In order to formulate the relationship between these two variables, $\frac{\partial B}{\partial u}$ at $O$ is calculated by numerical derivation, and the results are given in Fig. 5 as a function of $\epsilon$ . It should be emphasized here that because the magnetic propulsion at the micro-scale is considered in this study, the center point $O$ was taken as the point of interest since it is considered as the position of the microrobot whose dimensions are assumed to be as small as a point compared to the size of the electromagnets. The representative mathematical correlation is found by using exponential curve fitting as in Eq. (6) ( $R^2=0.99$ ).

(6) \begin{equation} G_o(\epsilon )=18.2(1-e^{-0.0368\epsilon }) \end{equation}

where $G_o(\epsilon )=\frac{dB(\epsilon )}{d_u}$ is the function of the normalized magnetic field gradient $\frac{\partial B}{\partial u}$ at the center point $O$ along $u$ axis based on $B_o$ (5.8 mT in Fig. 4). The volume $V$ and the magnetization vector $\vec{M}$ can be expressed with a structural constant $C_\mu =V|\vec{M}|$ since they are characteristic properties of the permanent magnet of the microrobot. Consequently, the force equation can be rewritten in a simplified form as in Eq. (7).

(7) \begin{equation} \vec{F}=C_{\mu } G_o(\epsilon ) \end{equation}

Equation (7) and Fig. 5 indicate that the magnitude of the magnetic force $\vec{F}$ is proportional to the magnetic field gradient $\frac{\partial B}{\partial u}$ which can be controlled by changing the asymmetry parameter $\epsilon$ . Hence, the asymmetry parameter $\epsilon$ can be defined as a control variable of the magnetic force magnitude.

4.1. Control strategy and task space control

The proposed strategy for the microrobot position control is to apply a programmed sequences of an oriented 3D force along a desired route. The initial requirement for this strategy is to locate the microrobot on the center point $O$ of the symmetric configuration when no current flows through the coils. When the control system is activated (so the electromagnets), the error between the current and the following position on the sequence is calculated by using feedback from the task space vision system. The system then generates a 3D magnetic force output. As the force is defined by its 3D direction and magnitude, the end effector orientation (i.e., coaxis) and the electromagnet asymmetric position $\epsilon$ are updated accordingly. The microrobot then moves and loses its initial position in the coordinate frame of the end effector toward the following position on the sequence. The robot arm moves continuously to compensate the deviation of the microrobot position in that coordinate frame. This compensation is repeated with coaxis orientation at 5 DOF until fulfilling the designed path.

According to this strategy, the control problem of the proposed EMA system can be considered as a combination of a prime task and an auxiliary task, as shown in Fig. 6. The prime task objective is to compensate the position of the microrobot in the 3D workspace according to the programmed path, while the auxiliary task is to keep its position constant in the coordinate frame of the end effector (i.e, relative position to the electromagnets). Losing the defined relative position will shift the microrobot to another gradient region-namely, different magnitude of the magnetic force – which will cause losing control and sudden attraction of the microrobot to one of the electromagnets. To avoid that, the control system takes the actual position of the micro robot from the 3D vision setup and calculates the positioning error in the coordinate frame of the end effector. Afterward, the robot arm moves to compensate the relative positioning error. The speed of the relative compensation should be much faster than the motion of the microrobot to sustain the stability of the EMA system. The difference in operation speeds of the utilized robot arm controller (shown in a following section) guarantees that the microrobot positioning error with respect to the end effector is compensated approximately 20 times faster than microrobot positioning error in the task space.

Figure 6.

Block diagram of the control system of the untethered EMA system.

Furthermore, proportional task space position controller is implemented to operate with the given position reference and the microrobot position feedback provided by the vision system. The positioning error is calculated in the 3D space as a vector with magnitude and orientation. The error magnitude is taken as a reference for the magnetic force magnitude controller. The magnetic force magnitude controller implements proportional control to calculate the required force magnitude by using the experimental relationship between $\nabla \vec{B}$ and the asymmetry parameter $\epsilon$ . The orientation error, on the other hand, is provided to the magnetic force orientation controller. The orientation controller updates the three DOF rotational references of the robot manipulator’s end effector. The combined effect of the magnetic force’s magnitude and orientation controller induces the task space force to compensate the microrobot positioning error according to the reference path. It is worthwhile to mention that controlling multiple microrobots within the ROI at the same time with the same strategy is possible if their magnetic fields are not interfering with each other. However, simultaneous independent multitask microrobot control is not achievable with the current EMA configuration.

4.2. Low-level controllers

As described above, the control of the electromagnet coaxial positions is necessary for setting the magnitude of the magnetic force. The translational motion of the electromagnets is provided by a timer belt mechanism that is driven by smart servo motors. Due to the usage of the servo motors, the mentioned control problem is solved by parameter tuning of an internal PI loop of the actuators. Furthermore, the real-time magnetic field control in this study does not involve any dynamical current updating. Therefore, the objective of the current control is to eliminate the steady-state error in the presence of disturbances. A custom PI controlled current driver is designed and implemented for this purpose [Reference Akcura, Alasli, Cetin and Tamer36].

4.3. Task space vision feedback system

The feedback for the control is provided by a calibrated orthogonal stereo camera setup. The orthogonal setup consists of a static bottom view camera fixed under the test vessel and a side view camera fixed to the end effector (like an eye in hand configuration). The stereo reconstruction is carried out by taking homogeneous coordinate transformations that relate three frames: the end effector camera, bottom camera, and the side view camera, as illustrated in Fig. 7. The required input of the 3D position reconstruction is the position of the object in the image plane of both cameras. Consequently, the position of the microrobot should be tracked in an image flow during the actuation. The “boosting” algorithm is used for this purpose due to its applicability in real time tasks [Reference Grabner, Leistner and Bischof37]. The algorithm provides the position of a predefined rectangle in the image plane. Then, image segmentation routines and image moment calculations are carried out to define the necessary data of the image plane for the 3D reconstruction.

Figure 7.

The used coordinate frames of the task space vision feedback system in stereo reconstruction.

5.1. Actuator prototype

The actuator prototype consists of a robotic arm, structural base of the EMA system, two coaxial electromagnets, electromagnet motion controller, electromagnet current controller, and task space feedback vision system, as shown in Fig. 8. Kuka Agilus KR6 R900 manipulator with positioning accuracy of 30 $\mu$ m is used as the robotic arm. The structural base is a 3D printed bridge-like assembly made of ABS enforced with carbon fibers. The base holds the pair of electromagnets on two coaxial sliders and brass guide lines that align and allow the linear motion of them on the same axis. The electromagnets were fabricated according to an optimized design for the purpose of the proposed actuator. The design took into consideration the size constraints of the EMA system, which is set to 500 mm as the maximum total coaxial length limit including the total length of the electromagnets (300 mm) plus the maximum length of the ROI (200 mm), as discussed in details in Ref. [Reference Alasli, Çetin, Akçura, Kahveci, Can and Tamer35]. The electromagnet motion controllers control the position of each electromagnet with a motorized timer belt driven by Dynamixel XM430 Smart Servo motor. The servo motors have built-in functions for implementing different control routines; the PI controller mode is selected for this application. With this setup, the distance between the electromagnets can be arranged from 80 to 200 mm (D1 and D2: 40–100 mm). Furthermore, a positioning accuracy of 0.1 mm is achieved, and the system has a frequency of 16 Hz. The electromagnet current controllers are designed to eliminate steady state errors in presence of disturbances since the real time magnetic field control in this study does not involve any dynamical current updating. A PI controller is designed and implemented for this purpose with an observed bandwidth of 500 Hz. Two Basler cameras (with 16 mm lens and 1280 × 1800 image resolution at 25 Hz) providing planar position feedback with a resolution of 20 $\mu$ m (at a fixed distance of 120 mm) is integrated to the system as the key-components of the task space feedback vision system. The proposed algorithm of the task space motion controller is coded on a single board computer NVIDIA Jetson TK1 and operated at 20 Hz.

5.2. Electromagnetic field experiments

To verify the capability of the actuator in generating homogeneous/gradient in the magnetic field at the symmetric/asymmetric configuration of electromagnets, the same conditions of the numerical analysis were applied. The magnitude of the magnetic field along the central line was measured by AFW Bell 5180 Gausmeter equipped with STD18-0404 probe at equal distances of $D_1=D_2$ 50, 60, and 70 mm for the symmetric configuration at constant excitation current of 2 A. Afterward, the distance $D_1$ was kept constant at 60 mm, and $D_2$ was shifted to 70, 80, 90, and 100 mm to observe the magnetic field gradient at the asymmetric configuration. The measured magnitude of the magnetic fields on the central line corresponding to the two configurations is plotted in Fig. 9.

Figure 9.

Comparison between the experimental and numerical results of the magnetic field at (a) Symmetric configuration and (b) Asymmetric configuration.

Figure 10.

Magnetic field gradient at center point $G_u (\epsilon )$ with respect to asymmetry parameter $\epsilon$ (Experimental Results).

The experimental and numerical results around the center point $O$ (region of interest) are consistent, as can be observed from Fig. 9 for both the symmetric and asymmetric configurations. However, deviation increment is observed as the measurement point approaches the electromagnet. The deviation can be attributed to the losses caused by environmental conditions and the nonlinear effects arising from the core structure. Furthermore, the gradient at point $O$ can be calculated by using the results of Fig. 9(b). For that, a function $B(u)$ of the magnetic field magnitude on the central line was obtained as a polynomial regression at the second order for each $D_2$ . The gradient at point $O$ is then evaluated by taking the derivative of the polynomial at $O$ . The calculated values are given in Fig. 10 with respect to the asymmetry parameter ( $\epsilon$ ), which is 0, 10, 20, 30, and 40 mm for $D_2$ 60, 70, 80, 90, and 100 mm, respectively. The representative mathematical relationship is found by employing exponential curve fitting as in Eq. (8) ( $R^2=0.97$ ).

(8) \begin{equation} G_u (\epsilon )=14.1\left(1-e^{(-0.366\epsilon )}\right) \end{equation}

By observing the fitted exponential curve in Fig. 10 and Eq. (8), it can be seen that it is in the form $K(1-e^{-\kappa \epsilon })$ , which is characterized by two parameters $\kappa$ and $K$ . The parameter $\kappa$ defines the spatial dependency of the function output on the asymmetry parameter $\epsilon$ , while K acts as a scaling factor. By comparing Eqs. (8) with (6) (of the numerical results), it can be noticed that $\kappa$ in both equations are close to each other (0.0366 in Eq. (6) and 0.0368 in Eq. (8)). The similarity in $\kappa$ indicates that the spatial dependency of the magnetic field gradient has same characteristics in both numerical and experimental studies. Furthermore, the scaling factors can be better compared when they are normalizing based on $B_o$ (5.8 and 5.6 mT for the numerical analysis [Fig. 4] and the experimental measurements [Fig. 9] values, respectively). The normalized $K$ in Eq. (8) (14.1) is smaller than the one in Eq. (6) (18.2) due to the weakening of the gradients with distance to pole’s face and the real world implementation.

5.3. Magnetic force experiments

A cylindrical magnet with a diameter of 1 mm and height of 1 mm is used as a microrobot prototype for this set of experiments. The remanence of the magnet is 1.28 mT, and the calculated magnetic dipole ( $C\mu$ ) value is $8.00\times 10^{-7} A.m^2$ . The microrobot motion is observed with the cameras, while it is moving in a fluid with a viscosity of 557 mPa.s inside the test vessel.

5.3.1. Constant force (gradient) experiments

As observed from Eq. (6), the relation between the magnetic field gradient $\nabla _u\vec{B}$ and the force acting on the magnetized part of the microrobot is proportional. Based on this, a set of experiments to analyze the microrobot motion under the effect of constant force were performed. The asymmetry parameter ( $\epsilon$ ) is changed from 10 to 40 mm with 10 mm steps at initial $D$ = 60 mm, while the relative distance of the microrobot to the electromagnets is kept constant. The experiments are carried out at open loop control and the force was applied for 100 s. The results of tracking the position of the microrobot with time at different $\epsilon$ are shown in Fig. 11. As it can be noticed from the slopes of the lines, the microrobot reaches a constant velocity at the steady state. Furthermore, increasing the gradient ( $\epsilon$ ) results in an increment in the steady-state velocity.

Figure 11.

Motion of the microrobot under the effect of a constant force.

Physically, microrobot velocity reaches the steady state when the electromagnetic force is at equilibrium with the drag force by compensating the gravitational force with buoyancy force and neglecting the lift force since the microrobot is considered perfectly symmetric. By exploiting this principle, the capability of the proposed actuator of manipulating microrobot at constant force can be evaluated. The steady state velocity and the Reynolds number for each experiment are first calculated and given in Table I. Since the Reynolds numbers are smaller than 1, drag force model (given in Eq. (9)) is applied according to Stokes’ law.

(9) \begin{equation} F_c=C_f.v \end{equation}

where $C_f$ is the drag coefficient that depends on the geometry of the microrobot and the viscosity of the fluid and $v$ is the measured velocity of the microrobot. By applying Newton’s second law, the equation of the motion can be written as in Eq. (10).

(10) \begin{equation} m\dot{v}=C_\mu G_u(\epsilon )-C_f.v \end{equation}

At steady state, the time derivative term is zero and a proportional relationship between the velocity and the magnetic gradient is obtained as a constant ratio $\frac{C_\mu }{C_f}$ , which is defined only by the characteristics of the fluid and microrobot. Since the experiments are performed with the same fluid and same microrobot, the ratio $\frac{C_\mu }{C_f}$ should be constant. Consequently, the actuator capability of controlling the force magnitude can be evaluated by calculating this ratio for each $\epsilon$ , as given in Table I. Based on the experimental results, the average of the ratios is 380.8 $\frac{mT.s}{mm^2}$ with standard deviation of 7.2 $\frac{mT.s}{mm^2}$ . Since the standard deviation corresponds to 2% of the average, the ratios are assumed to be constant, and it can be said that it is possible to sustain constant magnitude of the induced force during actuation via the proposed actuator.

Table I.

Experimental results and calculated parameters of the constant force experiments.

5.3.2. Transient force (gradient) experiments

The capability of generating transient force is tested with sinusoidal force (gradient). To apply a sinusoidal gradient on the microrobot at $O$ of the form $G_u(t)=|G_u(t)|sin(2\pi ft)$ , the asymmetry parameter ( $\epsilon$ ) should be given as a generic sine wave input to the motion controller of the electromagnets. For that, the time dependent form of $\epsilon$ , Eq. (11), was first obtained from the inverse of the experimental-based $G_u(t)$ function (Eq. (7)). $\epsilon (t)$ data according to the frequency was then generated and fed to the two smart servo motors who drive the coaxial electromagnets. The applied frequency was changed from 0.01 to 0.125 Hz, while its amplitude is kept constant at 0.025 T/m. The duration of each experiment is set to 6 times the period of the input. The last three oscillations were taken to ensure analysis at the steady state.

(11) \begin{equation} \epsilon (t)=-\frac{1}{0.366}\frac{log(1-G_u (t))}{14.1} \end{equation}

Figure 12 presents the microrobot responsive motion, the actuator commands to the left and right motors of the electromagnets, and the motion of the end effector of the robotic manipulator at a 0.01 Hz sinusoidal input. As expected, the microrobot makes steady state simple harmonic motion on the central line under the effect of the sinusoidal force. Furthermore, the synchronized motion of the end effector for compensating the microrobot motion in order to keep it at the center point (i.e., constant $\nabla \vec{B}$ ) can be also realized.

Figure 12.

Steady state gain with respect to frequency.

Table II summarizes the obtained results of the sinusoidal magnetic gradient (force) experiments. The results present the frequency dependent change in the steady state gains. The gain is calculated as the ratio of the sinusoidal steady state microrobot velocity magnitude (output) to the sinusoidal gradient magnitude (input). The gains can be considered as a representative value of the mechanical admittance ( $v/F$ ) for the interaction between the actuator and microrobot since the magnetic force is proportional to the magnetic field gradient. The results also reveals that the bandwidth of the experimental actuator is 0.1 Hz.

5.4. 3D motion demonstration

Manipulation of a microrobot model on a helical route is demonstrated in Fig. 13 and the supplemental video. At initial stage, the programmed path is uploaded to the space task controller. The microrobot model is placed at the central point with respect to the robot arm and both cameras are activated to log the position of microrobot as the central position. Then, the electromagnetic actuator is activated and a magnetic force is applied on the microrobot. At the proceeding time instant, the control system receives the position of the microrobot while its moving from the cameras and calculates its displacement error with respect to the central position on the image plane. According to the programmed path, the space task controller applies proportional control to give an incremental position and direction command to the robot arm by scaling the observed motion in the image plane to the real world coordinates in order to compensate the displacement error. The last loop is repeated until the microrobot fulfils the whole path. Figure 13 presents the tracks of the microrobot, robot arm on XY plane (bottom camera), YZ plane (frontal camera) and the stereo reconstruction 3D image by the space task vision system according to the reference path. Moreover, the average of the absolute error according to the nearest point on the reference path and the standard deviation $\sigma$ are 190 and 101 $\mu$ m, respectively, as it can also be realized from the graphs. It is worth mentioning here that a part of the error is related with the robotic arm and image noise (see the supplemental video).

Table II.

Frequency response of the actuator (Magnetic gradient vs microrobot velocity).

Figure 13.

The tracks of the microrobot and the robot arm on (a) XY plane (bottom camera), (b) YZ plane (frontal camera), (c) the stereo reconstructed 3D path and (e) the absolute error according to the nearest point on the helical reference path.

The performance of the proposed system is comparable with other EMA systems in terms of positioning accuracy and velocity (Table I). For instance, the stationary EMA system, developed by Khalil et al., propelled a sperm-cell shaped microrobot at low magnetic field strength of 1.39 mT at an average velocity of 25 $\pm$ 10 $\mu$ m/s with an accuracy of 90 $\pm$ 40 $\mu$ m [Reference Khalil, Magdanz and Sanchez38]. Similar stationary systems also reported positioning errors of 100 $\mu$ m at 122 $\mu$ m/s for microsphere robot (r = 350 $\mu$ m) [Reference Ongaro, Pane, Scheggi and Misra3] and 310 $\mu$ m for 400–4000 $\mu$ m multi-microrobots [Reference Diller, Giltinan and Sitti39]. On the other hand, mobile EMA systems such as the ARMM system of Heunis et al. [Reference Heunis, Wotte, Sikorski, Furtado and Misra40], positioned a magnetic endovascular catheter at an accuracy of 2.09 $\pm$ 0.49 mm. Moreover, it is worth mentioning here the conclusion of Franz et al. of their comprehensive review of electromagnetic tracking systems in medicine. They stated that electromagnetic steering systems have an average accuracy of 1 mm in minimal interference environments. Nevertheless, any source of distortion can cause drastic accuracy drops [Reference Franz, Haidegger, Birkfellner, Cleary, Peters and Maier-Hein41].