Development of a novel nonlinear model and control strategy for soft continuum robots featuring hard magnetoactive elastomers

Magnetoactive soft continuum robots (MSCRs), capable of controllable steering and navigation, hold substantial promise for healthcare applications. However, advancements in MSCRs have been hindered by a limited understanding of MSCR dynamics and a lack of effective control methods. Addressing these gaps, this study presents a novel, time-dependent, and computationally efficient analytical model of MSCR, alongside a new optimal closed-loop control strategy for precise high-frequency trajectory tracking. A finite element (FE) model of the MSCR is initially developed, with its validity confirmed through rigorous laboratory measurements. Using the formulated FE model, a new and computationally efficient analytical model is subsequently developed to accurately predict the highly nonlinear response of MSCR. This model operates as a system of switched linear models, each of which is a reduced-order version of its corresponding high-order linear model extracted from the FE analysis. This innovative approach not only maintains the predictive accuracy of the FE model but also significantly reduces computational demands, operating in just a few seconds. The results highlight that the developed model can accurately predict the dynamic responses of the MSCR while significantly reducing the computational load by almost 80 orders of magnitude compared with the FE model on the same simulation platform. The proposed model has been effectively utilized to develop a novel optimal control strategy using the feedforward interval type-2 fractional-order fuzzy-PID method. A hardware-in-the-loop experimental test has been finally designed to demonstrate the superior performance of the MSCR under the proposed controller.


Introduction
In recent years, designs of soft continuum robots (SCRs) have emerged to achieve enhanced steering and navigation capabilities in constrained, confined and complex environments [1][2][3].These offer promising potentials for minimally invasive therapies (figure 1(a)), such as endovascular treatment of cardiovascular diseases.A number of SCR designs have evolved during the past decade particularly for a range of minimally invasive cardiovascular therapies [4][5][6].In such treatments, a thin tube or a passive guidewire, known as a catheter, is manually guided by radioscopic imaging into the heart chamber up to the targeted locations, e.g. in front of blood vessel branches or lesions tissues, to ablate the lesion [2].Although the commercially developed SCRs have proven to be effective in reducing the associated risks and the side effects compared to the conventional methods, the current designs are known to impose notable operational challenges.These are mostly associated with their actuation mechanisms, such as pulling mechanical wires, inflating pneumatic or hydraulic chambers, and embedding rigid magnets for manipulation, apart from difficulties in accessing small branches [7].These, together with relatively longer operation times, make the current designs unfavorable for numerous biomedical tasks [1].
In recent years, magnetoactive soft continuum robots (MSCRs) with multimodal locomotion designs have been proposed as alternatives for the distal portions of the SCRs, while circumventing some of their performance limitations [7][8][9][10][11].This class of soft robots generally comprises magnetoactive elastomers, in which magnetized hard micron-sized particles such as neodymium-iron-boron (NdFeB) are impregnated into a soft elastomeric matrix such as silicone rubber or polydimethylsiloxane (PDMS) [12][13][14].Figure 1(b) illustrates an elastomeric cantilever MSCR with embedded NdFeB particles representing the distal portion of the SCR.The MSCRs offer the superior potential for wireless active navigation and steering in unstructured and confined environments under a relatively low-level magnetic field stimulus.Moreover, these soft robots exhibit superior bending flexibility allowing access to very small branches, as seen in (figure 1(c)) [1]. Figure 1(c) shows a MSCR consisting of a magnetoactive distal portion that undergoes large bending deformation under the applied magnetic field, while a non-magnetic segment controls the advancement and retraction of the MSCR.Moreover, using ferromagnetic microparticles as an actuation source can permit miniaturization of the MSCRs to sub-millimeter scales, which would be favorable for more precise, complex, and repeatable delicate biomedical tasks in vivo/vitro settings, especially for minimally invasive therapies [7,8,13,15].The magnetic body force and torque developed by the embedded hard microparticles under an external magnetic field also enable the MSCR to operate wirelessly within enclosed environments, e.g. the human body.
The rapid advances in 3D and 4D printing technologies [16,17] together with developments in novel magnetic circuits and solutions [18][19][20] have contributed to developments in the fabrication process and actuation mechanisms for MSCRs [19].The formulations of effective dynamical models and practical closed-loop control algorithms, however, are still in their infancy, which are vital for developments and real-life practical applications of MSCR for minimally invasive surgical tasks.Only limited efforts have been reported for dynamic modeling and closed-loop control of MSCR [21][22][23][24][25]. Zhao et al [24] proposed an analytical model employing linear theory to estimate the bending angle of an MSCR in the static state.The MSCR considered in the study comprised a soft beam with hard-magnetic particles, while the developed analytical model was limited only to a small deformation regime.The study also developed a finite element (FE) model of the MSCR beam using ABAQUS considering large magnitude bending deformations, the validation of which was examined via laboratory measurements.Chen et al [25,26] formulated a mathematical model for estimating the static deformation of a hard-magnetic soft beam considering large deformations.The predicted responses revealed good agreements with experimental data under high magnetic field intensities, while considerable deviations between the model-predicted and the measured responses were noted under lower magnitudes of the magnetic field.Dadgar-Rad and Hossain [27] conducted a time-dependent finite deformation analysis of a soft beam with hard-magnetic particles subject to a magnetic loading.In a recent study, Kadapa and Hossain [28] developed a constitutive model of an iron-filled magneto-active polymer using a unified numerical approach.
The modeling efforts reported in the aforementioned studies have been either limited to static analyses or involve computationally demanding numerical analyses.Such models are not thus suited for developments in control algorithms, which are vital for realizing closed-loop navigation and tracking performance.The tracking performance characteristics of soft magnetic robots have been mostly studied using open-loop control strategies or simple closed-loop methods, such as P, PI, and proportional-integral-derivative (PID) [22,23,29].Apart from their poor accuracy, such methods yield poor tracking performance in the presence of disturbances and uncertainties.Kim et al [18] developed an open-loop strategy via telerobotic for minimally invasive therapies to steer a magnetic SCR into narrow pathways of neurovascular phantoms.Zhang et al [30] developed a small-scale SCR with a largeangle steering capacity.In order to steer the robot into desired sites, an open-loop control strategy was formulated to adjust the external magnetic field.Liu et al [31] developed a remotely controlled MSCR comprising NdFeB particles dispersed in a PDMS matrix to steer the robot in a bifurcated environment.The study employed a number of predefined magnetic field stimuli of different magnitudes and directions to control the steering of the MSCR.
Motivated by the limitations of conventional SCRs, design challenges, and superior potential of the MSCRs, this study presents a computationally efficient time-dependent model of an MSCR together with a high-performance closed-loop control algorithm to achieve accurate tracking of predefined trajectories even at relatively high-frequencies.A high-fidelity FE model of an MSCR, fabricated in the laboratory, is initially developed and analyzed under a perpendicular uniform magnetic field using COMSOL Multiphysics ® software.The validity of the proposed FE simulation is investigated using the measured deformation responses of the MSCR under different magnetic excitations.Subsequently, using the developed FE model, a computationally efficient time-dependent nonlinear model is formulated by combing a number of weighted linear models, each of which is a reduced-order version of its corresponding high-order linear model extracted from the FE analysis.The validity of the reduced model is examined by comparing predicted deformation responses with those obtained from the FE model.A novel feedforward interval type-2 fractional-order fuzzy-PID (FFIT2FOF-PID) control algorithm is finally designed and implemented to achieve improved tracking performance and robustness against uncertainties by taking advantage of an interval type-2 fuzzy logic (IT2FL) network and fractional operators.The parameters of the developed control algorithm are optimally tuned using the general modified Cuckoo optimization algorithm (GMCOA) considering a set of predefined trajectories.The effectiveness of the proposed algorithm is demonstrated in terms of steering and tracking performance of the MSCR.The developed model and control strategy may provide essential guidance for the future development and practical implementation of this emerging technology in real-life applications.The remaining manuscript is organized as follows: The governing equations of the MSCR are initially presented in section 2. Section 3 describes the formulation of the MSCR using the weighted reduced-order linear models extracted from the developed FE model.The proposed control strategy is subsequently formulated in section 4, while the materials and methods are presented in section 5. Simulation and experimental results are finally discussed in section 6, and the conclusion of the current study is provided in the last section.

Model formulation
The MSCR, shown in figure 1(b), is made of a non-magnetic silicone-based elastomeric substrate impregnated with hard micron-sized magnetic particles (NdFeB).This MSCR is modeled as a magneto-rheological (MR) elastomeric cantilever beam subject to a uniform magnetic field normal to the robot's length.The residual magnetic flux density induced in MSCR is assumed to be aligned along the x-axis in the static equilibrium configuration.
The magneto-sensitive composite structure of the MSCR together with the applied magnetic field constitutes a coupled magneto-mechanical problem.The governing equations of motion of the coupled system can be described by combining elasticity and Maxwell equations in the presence of an external magnetic load [32]: where E, F B , I D , B, H, ρ s , u s and σ s are the electric field, body force, current density, magnetic flux density, magnetic field intensity, mass density, displacement of the robot structure, and Cauchy stress tensor of the solid medium, respectively.The magnetic flux density 'B' and the current density 'I D ' are related to magnetic field intensity 'H' and electric field 'E', respectively, considering the following constitutive relationships: where µ 0 denotes the vacuum permeability, M and σ are the remanent magnetization and conductivity of the medium, respectively.The elastomeric structure is modeled using the widely used Neo-Hookean hyperelastic material model to describe its mechanical behavior.Assuming incompressibility of the elastomeric material of the MSCR's substrate, the strain energy function 'W' for incompressible Neo-Hookean material in the three-dimensional space can be expressed as [33]: where C 01 and I 1 are the material constant and the first strain invariant or trace of the left Cauchy-Green deformation tensor 'G', respectively.In a continuum medium, strain invariants (I 1 , I 2 , and I 3 ) can be obtained from the left Cauchy-Green deformation tensor 'G', which is directly related to the finite deformation gradient tensor 'F', such that [34]: In the above relations, J denotes the Jacobian matrix or the determinant of the finite deformation gradient tensor 'F' and λ 1 , λ 2 , and λ 3 are, respectively, the stretch ratios in the length, width, and thickness directions, which are evaluated from: where w, l, and h are the length, width, and height of the rubber-like substrate in which notations with the subscript ' 0 ′ denote their initial state.For the incompressible material subjected to uniaxial tension, the stretch ratios are related as [33]: Equations ( 7) and ( 11) yield following expression for strain energy density as a function of λ 1 : Considering the Cauchy stress for the incompressible Neo-Hookean hyperelastic material model and uniaxial stretch, the stress in the pulling direction is subsequently obtained as: where k is the shear modulus or the second lame parameter.In this study, the material constant 'C 01 ' is identified from the measured uniaxial loading test data, which is subsequently used to estimate the elastic modulus of the elastomer 'E m ' as: where ν p denotes the Poisson's ratio of the incompressible elastomeric substrate.
The magneto-elastic and magnetic behaviors of the MSCR are coupled by the Maxwell surface stress developed at the boundary between the magnetic and elastic domains, can be described by the following relations [32]: where n is outwards normal to the surface bounding the magnetic and non-magnetic domains, and σ 1 and σ 2 denote the stress tensors at the interface attributed to magneto-elastic behavior of the MSCR body (subscript 1) and the magnetic domain (subscript 2), respectively.The subscript ' H ' in the above relations denotes the full Maxwell stress tensor on the boundary of the MSCR body due to the applied magnetic field, given by [32]: where δ ij is Kronecker's delta.The above formulations can be solved to determine the deformation response of the MSCR body under the application of an external magnetic field.In this study, a FE model of the MSCR is developed in the COMSOL Multiphysics platform to study the magneto-elastic responses of the MSCR subject to a magnetic flux density 'B', as described in section 6.The flux density input also relates to the electromagnet coil current 'I', which can also be related to voltage output of the controller considering a simple R-L circuit for the coil.

Reduced-order linearized model
The magneto-elastic analysis of the MSCR is a complex nonlinear dynamic problem.A magneto-elastic FE model can be effectively formulated to accurately predict the response behavior of MSCR under varied applied magnetic fields.Such a high-fidelity model is, however, computationally very heavy and thus not suitable for the design of the controller.In this study a novel reduced-order model of the MSCR is derived.
The developed model practically enables the design of the controller for enhancing the tracking performance using the trajectory piecewise linear (TPWL) approach [35].The TPWL technique seeks to estimate the nonlinear magneto-elastic behavior of the MSCR over a broad range of input, considering a combination of weighted linearized models relying on multiple linearization points evaluated using the full FE model.The nonlinear dynamic responses of the MSCR may be expressed in the general form, as: where xϵR n and yϵR n describe the states and outputs, respectively, and uϵR p denotes the control inputs.Moreover, f : R n → R n and g : R n → R q are nonlinear functions.Assuming negligible contributions due to high-order terms, the linearized model corresponding to a given state (x (i ) , u (i ) ) can be expressed as: where, A i and C i , respectively, denote the Jacobian matrices of f (•) and g (•) estimated at (x (i ) , u (i ) ), and B i = ∂f ∂u evaluated at (x (i ) , u (i ) ).Assuming the steady-state condition, f x (i ) , u (i ) is considered to be zero in the above relations.Letting x = x − x (i ) , ũ = u − u (i ) and ỹ = y − g x (i ) , the linearized model can be expressed as: The above linear model is considered valid in the vicinity of the chosen operating state x (i ) .The effectiveness of the linearized model can be, however, enhanced by considering a set of linear equivalent models corresponding to different operating points along a defined trajectory.This will permit an accurate estimation of the nonlinear system's response over a broad range of inputs.Since the current state 'x' of the designed linear model lies near the defined trajectory, the closest linear model (A i , B i , C i ) approximately mimics the nonlinear response around the current state.The distance d i between the current state of the linearized model and the linearization point is estimated by the Euclidean norm, x − x (i )  2 .The response of the nonlinear system is then captured through appropriate switching among different linearized models.
Consider a two-dimensional (2D) trained trajectory, shown in figure 2, with various linearization points x (i ) .A weighted linearized model is formulated considering a combination of linear models corresponding to different linearization points x (i ) , u (i ) ; i = 1, 2, . . ., m , where m denotes the number of linear models or operating points considered.The weighted equivalent model can thus provide an accurate estimate of the nonlinear responses, and is given by: where, ω i (i = 1, 2, . . ., m) are the normalized weighting factors ( ) associated with the linearized models corresponding to points (x (i ) , u (i ) ).The weighting factors are then identified at each time step by evaluating the distance error between the current state and linearization points ' 2 ', using the following algorithm, leading to the identification of the dominant linearized model associated with the linearization points x (i ) at each time step [35].
In the algorithm, β is a positive control parameter, which affects the selection rate of dominant linear models.As it can be realized, the linear model with higher d i has a lower weighting factor, thus less dominant.
The linear model derived on the basis of the results obtained from the FE model generally comprises a large number of degrees-of-freedom (DOF).The modelbased simulations, and the controller design process thus becomes computationally demanding, and generally infeasible for real-time tracking control.Development of a reducedorder model is highly desirable not only for reducing the computational demand but also for facilitating the controller design.It is, however, vital that the reducedorder model (ROM) yields the system's responses with reasonable accuracy.In this study, the modal truncation algorithm, sssMOR [36], is employed to reduce the order of the extracted linear system.Considering equation ( 19), the reduced-order linear model may be expressed in the following form: where zϵR v denote the states of the reduced model, A r ϵR v×v and B r ϵR v×r .Here, a transformation matrix 'V' exists that maps the states of the original linear model to those of the ROM.Using equations ( 20) and ( 21), the ROM is formulated as follows, where w ∈ R r denote the outputs: Owing to its relatively lower dimensionality, the ROM can provide accurate predictions of the nonlinear system's responses in a computationally efficient manner.Moreover, the reduced-order model can facilitate the development of advanced control strategies.

Feedforward interval type-2 fractional-order fuzzy-PID controller
A high-performance closed-loop control strategy based on IT2FL is designed.The objective is to minimize the tracking error of the proposed MSCR in real-time.The tracking error 'e' is defined as the deviation of the tip angle 'θ' of the MSCR with respect to the desired angle 'θ d '.The tracking error and its fractional time derivative 'D α e' can be expressed as: Remark 1.In this study, the Grunwald-Letnikov definition is used for fractional-order derivative and integral operators [37].The proposed IT2FL network is designed considering two inputs (x f1 , x f2 ) based on the tracking error and its fractional time derivative as: where K P and K D are constant coefficients, which are determined using the GMCOA [38,39].The inputs x f1 and x f2 are partitioned by N type-2 fuzzy sets, labeled as, Fi 1 (i = 1, 2, . . ., N) and Fj 2 (j = 1, 2, . . ., N), respectively.Since there exists a rule for each possible combination of N sets of the two input domains (x f1 and x f2 ), the rule-base network comprises M = N 2 rules.The kth k = N (i − 1) + j, k ∈ 1 M rule 'R k ' of the IT2FL network is assumed to be of the form: where k zk ] denotes the consequent membership functions (MFs), and z − k and zk are crisp numbers that are located at a distance ν away from either side of the center of symmetry.The firing set associated with the kth rule, F k , is defined considering the interval set given in equation (28).Each input in the aforementioned rules takes seven IT2FL MFs, including two Sigmoid and five Gaussian type functions, which are defined in the following relations: where, In the above relations, c A and c B are the MFs centers associated with the first and second inputs of the designed IT2FL network, respectively.σ uA , σ lA , σ uB , σ lB , η lA , and η lB are the coefficients of the Gaussian and Sigmoid functions, whose values are determined using the optimization process.µ F1 and µ F2 denote the membership functions of the transverse angle error and its fractional rate, respectively, in which μF1 and μF2 are upper MFs, and µ F1 and µ F2 are lower MFs.Since returned values are type-2 fuzzy, the type of returned membership values is reduced to type-1 fuzzy using the center-of-sets typereducer method defined in equation (39) [40,41], such that: where the defuzzified outputs y l and y r are obtained from:  Error rate (x f2 ) In the above relations, R and L are the switching points, which are estimated via the iterative Karnik-Mendel (KM) algorithm [42].The mean of the outputs y l and y r is subsequently taken as the crisp output of the IT2FL network.The structure of the designed IT2FL network is depicted in figure 3, while the corresponding rules are summarized in table 1.
A novel FFIT2FOF-PID control algorithm is subsequently designed, as shown in figure 4, to achieve improved tracking performance and robustness against uncertainties.The algorithm employs velocity and acceleration feedforward compensations to further mitigate tracking errors [43].The overall control law U (t) of the FFIT2FOF-PID control algorithm is expressed as: where y (t), the IT2FL component of the control algorithm, is a nonlinear function of the error and its fractional-derivative, given by: where δ, K 1 and K 2 are constant parameters, which are estimated using the GMCOA.

Calibration of control algorithm's parameters
The parameters of the designed FFIT2FOF-PID control algorithm are tuned using GMCOA with the aim of minimizing an objective function (OF), defined as a weighted sum of some well-known indicators.These include integral absolute error (IAE), integral time absolute error (ITAE), and control energy factor (CE).The OF, together with the considered indicators, are expressed below: where t and t e denote the current instant and the entire test time, respectively.The time variable t in the considered OF above minimizes the oscillation chance in the system output and effectively reduces the closed-loop system's settling time [44].The optimization function, defined in equation (44), is subject to a number of limit and equality constraints in order to limit the number of design variables and achieve rapid convergence of the solutions, as defined in equation (48).Moreover, the IT2FL membership functions, F 4 1 and F 4 2 , are placed at zero positions, while the remaining MFs are positioned on two sides of these two membership functions.

Sample preparation and characterizations
The MSCR was fabricated by dispersing hard magnetic particles NdFeB (MQFP-B+, Neo-Magnequench) with a diameter of 5 µm in a silicone-based matrix.The SE1700 (baseto-curing agent mass ratio 10:1, Dow Corning) was mixed with Ecoflex 0030-Part A (Smooth-On, Inc) with a mass ratio of 1:1 to achieve desirable modulus of elasticity, while the volume of iron particles was kept at 17%.It is important to mention that a volume fraction range of 15%-20% is commonly used in the literature for magnetoactive soft actuators to achieve large deformation at low-intensity magnetic fields [7].Initially, the NdFeB particles were manually mixed into the resin using a spatula.Subsequently, this mixture was thoroughly blended with the particles in a Thinky ARV-200 vacuum mixer.The process was conducted under a vacuum pressure of 27 in Hg (91.4 kPa) and at an angular speed of 2000 rpm for a duration of 3 min to ensure homogenous mixing.Given the high viscosity of the resin mixture, which presents a challenge in removing air bubbles, multiple short-duration vacuum mixing sessions were employed instead of a single, prolonged session.This approach was strategically chosen to effectively eliminate air pockets from the mixture, while also mitigating the risk of partial curing that could arise from excessive heat generated during a longduration mixing process.An exact rectangular mold, measuring 67 × 65 × 10 mm 3 and designed in SOLIDWORKS 2021 to accommodate ten MSCR samples, was subsequently fabricated using an SLA 3D printer (Formlabs Form 3+) with an adjustable accuracy of 25 µms.Following fabrication, this mold was impregnated with a release agent (Ease Release 200, Smooth-On Inc) to facilitate sample removal.The prepared mixture was subsequently poured into the mold and permitted to cure in an oven at a temperature of 93 • C (200 • F) for about four hours.The cured sample, with desired dimensions of the MSCR (table 2), was removed from the mold using a razor cutter.Figure 5 presents the designed and fabricated mold and the fabricated MSCR sample.The prepared MSCR sample was subsequently magnetized in the direction of the robot's length using a strong magnetic field (1.6 T), provided by a variable airgap electromagnet (DXSBV-100, Dexing Magnet Tech.Co., Limited) equipped with 4 cm poles.The remanent magnetization of the sample was computed using the datasheets provided by Neo-Magnequench [45], which included the hysteresis demagnetization loop and the saturation level data.The fabrication process is illustrated in figure 6.It is noted that the distribution of iron particles was examined via a scanning electron microscope (Hitachi S-3400N), which is shown in figure 7(a).The micrograph reveals NdFeB particles' homogenous and random distribution within the matrix.Laboratory experiments were performed to identify mechanical property of the fabricated sample.The magnetized sample was subjected to tensile loading at a rate of 30 mm min −1 until failure via a uniaxial loading test system (F1505-IM, Mark-10).To ensure reliability, these tensile loading tests were conducted multiple times.The measured data were used to obtain stress-stretch characteristics of the sample, as shown in figure 7(b).The experimental data were subsequently fitted to the Neo-Hookean hyperelastic model, described in equation ( 13), to identify model parameters, as seen in figure 7(b).It is important to mention that for each tensile test, the modulus of elasticity was calculated as described above, and the value reported in the paper represents the median of all these measurements.The mechanical and magnetic properties of the MSCR sample are also summarized in table 2.

Measurements of magneto-elastic responses
An experiment was designed to characterize magneto-elastic responses of the MSCR subject to various magnetic field stimuli.The measured responses were used to examine the validity of the developed FE as well as ROM and also to assess the effectiveness of the proposed FFIT2FOF-PID control algorithm.The experiments were performed under a uniform magnetic field ranging from 0.5-32 mT, which was generated by a variable airgap electromagnet (DXSBV-100, Dexing Magnet Tech.Co., Limited) powered by a linear amplifier (LVC2016, AE Techron).The generated magnetic field around the MSCR was measured using a portable gauss meter (Parker 5180, F.W. BELL) and acquired in a data acquisition card (M Series-NI PCI-6251, National Instruments).The magnetoelastic response of the tip of the MSCR (angular deflection)  was captured by a high frame rate camera (Blackfly BFS-U3-17S7C-C, FLIR) together with the image processing algorithm programmed in Python language.A hardware-in-the-loop (HiL) was subsequently designed to measure the real-time actuation of the MSCR with the proposed FFIT2FOF-PID control algorithm implemented on the target computer using MATLAB-Simulink Real-Time (R2017a) with a sample time of 0.01 s.The experiment design for real-time actuation and control is illustrated in figure 8.

Results and discussions
The magneto-elastic responses of the MSCR, depicted in figure 1(b), are initially evaluated via a 2D FE model  developed in COMSOL Multiphysics (5.6) software.The model is formulated using the geometrical and material parameters presented in table 2. Finite deformation kinematics is implemented to account for large deformations using the Eulerian-Lagrangian formulation.The Neo-Hookean hyperelastic material model, described in section 2, is used, whose parameters were identified from the uniaxial tension test data (figure 7(b)).The MSCR is modeled within a sufficiently large airbox to ensure uniformity of the applied magnetic field surrounding the MSCR.The triangular elements of quadratic order are used to model the airbox, and the MSCR's body is modeled using the mapped mesh of quadratic order, which benefits from a uniform distribution pattern, as shown in figure 9(a).It should be noted that the considered mesh profile comprises 183 400 DOF.The validity of the chosen mesh profile was ensured through convergence tests under a uniform magnetic flux density of 15 mT until the reaction force at the lower fixed node boundary of the MSCR converged to less than five percent in subsequent simulations.
Figure 9(b) illustrates the steady-state deformation responses of the simulated MSCR under different magnetic field intensities ranging from 0.5 mT to 30 mT.Moreover, the experimentally measured deflection responses under different magnetic loadings (2 mT-25 mT) are shown in figure 10.
The results suggest nearly linear behavior of the MSCR at low magnetic field intensities and highly nonlinear behavior as well as saturation tendency under magnetic field intensity in excess of 15 mT.The validity of the model is examined by comparing the transverse angle of the tip of the MSCR with the measured data under different magnitudes of magnetic flux density, as seen in figure 11.The peak tip deflection responses obtained from the experiments and the model are also compared in table 3. Comparisons show reasonably good agreements between the FE model-predicted and measured responses.The peak error was observed as 2.7%, which occurred at very low magnetic field intensities, while the peak error under medium to high magnetic flux densities was less than 1%.

Reduced-order model simulations and verifications
The effectiveness of the developed reduced-order model of the MSCR is investigated by comparing its response behavior with that of the FE model for a defined trajectory.As mentioned before, the accuracy of the reduced model strongly relies on the number of operating linearized points (NOLP) or linear system models considered, as it is evidenced in equation (22).The model accuracy also depends on the order  of the reduced-order linear system (ORLS).It is, however, desirable to derive an effective ROM on the basis of the least number of linearization or operating points and least order.In this study, the reduced-order models are formulated considering 12, 18, and 36 operating points extracted from the FE model responses.Moreover, each reduced-order linear system model is designed with three different orders, namely, 2, 4, and 6.A total of 9 ROMs of the MSCR are thus formulated, which are labeled as (A.j, B.j, and C.j; j = 1, 2, 3) in table 4. The relative response prediction accuracy of the models is assessed considering a trajectory governed by 5 mT harmonic magnetic flux density stimulus at a frequency of 0.5 Hz.
Figures 12(a)-(c) compare tip deflection responses of the reduced-order models with those of the FE model.The comparisons suggest that the reduced-order models yield reasonably good estimations of the nonlinear responses, irrespective of NOLP and ORLS combinations considered.The deflection response of the reduced-order model with NOLP = 36 and ORLS = 6 exhibits very close agreement with that obtained from the FE model.The relative effectiveness of the reducedorder models is further assessed in terms of root mean square error (RMSE) and standard deviation (STD), which are reported in table 4. The results show that the ROM with NOLP = 12 and ORLS = 2 (A.1) yields RSME of 0.0240, which reduces to 0.0194 when ORLS is increased to 6 (A.3).The RMSE of the model with NOLP = 36 and ORLS = 6 (C.3) is observed to be substantially lower (0.0093).The STD of the error also shows a pattern similar to the RMSE.The computational efficiency of the reduced-order models is substantially superior to that of the FE model.The simulations of the ROMs were conducted using MATLAB-Simulink software (R2020b) on a dual-CPU computer system.Each Intel Xeon Silver 4210 R CPU features ten cores and twenty threads.The system is further equipped with 64 GB of RAM.The simulation of the FE model was also performed using COMSOL Multiphysics (5.6) installed on the same computer.Total simulation time was limited to 5 s for both the reduced-order and FE models.The run time associated with all the models, summarized in table 4, clearly shows a substantially lower computational demand for the linear models.For instance, for ROM with NOLP = 36 and ORLS = 6, the run time is nearly 6 s compared with 5100 s (85 min) required for the FE model.Owing to its ability to predict the MSCR responses accurately with superior computational efficiency, the ROM with NOLP = 36 and ORLS = 4 is used for the design and assessments of the developed FFIT2FOF-PID controller.

Calibration of the FFIT2FOF-PID controller
The parameters of the proposed FFIT2FOF-PID controller were optimally tuned using GMCOA due to its easy implementation, fast convergence, and capability to capture global optimum [38] in order to enhance the tracking performance of the MSCR with the reduced control effort.The developed ROM (case: C.2 in table 4) was effectively utilized to design the controller for a predefined trajectory.The control voltage obtained from the control algorithm (figure 4) serves as the input of the cascaded system considering the R-L model of the electromagnet.The resistance (R) and inductance (L) of the electromagnetic system were measured in the laboratory as 8.9 Ω and 3 H, respectively.
The parameters of GMCOA are summarized in table 5.These parameters, including the initial number of cuckoos (NICuckoo), the number of cuckoos in infinite iteration (NICuckoo inf ), the maximum number of allocated eggs to each cuckoo (Egg Max ), motion coefficient (MC), the impact factor of cuckoo fitness (C R ), the distribution angle (θ C ), and the positive constants λ C , Ω, and t C were determined by trial and error considering the best obtained solution in the least number of iterations with a high convergence rate.Due to the random  nature of GMCOA, the designed optimization problem was solved five times in an attempt to report the best pseudoglobal solution as FFIT2FOF-PID's parameters.The convergence graph of the GMCOA for the best solution and the normalized IT2FL surface rule after the optimization process are depicted in figures 13(a) and (b), respectively.The optimized IT2FL membership functions are also shown in figure 14.

Performance assessment of the FFIT2FOF-PID controller
The tracking performance of the optimally tuned FFIT2FOF-PID controller is assessed through HiL, as described in section 5.A total of 6 different target trajectories are considered for evaluating the control performance, which are given in equation (49).These include the harmonic trajectories in the 0.5-3 Hz frequency range with amplitude of an π 6 and a complex harmonic trajectory (combination of several harmonics).A trajectory with a relatively higher tip angle amplitude of π 4 is also considered to assess the control performance for dynamic trajectories of different amplitudes and frequencies.All cases are experimentally captured for a duration of 30 s using the Simulink Real-Time (R2017a).The proposed controller performance is also compared with those of the Fuzzy-PID (FPID) and PID methods, and its relative merits are discussed below.
Figure 15 illustrates the tracking performance of the proposed controller for the six desired trajectories defined in equation ( 49).The control effort in each case is also presented in figure 16 in terms of the required control voltage.Figures also show the tracking as well as control effort performances of the FPID and PID control methods for the sake of compassion.As demonstrated in figure 15, the tracking error of the FFIT2FOF-PID algorithm is considerably smaller than those of the FPID or PID methods for all of the considered cases (D.1 to D.6).The tracking performance of the FFIT2FOF-PID algorithm is also shown in supplementary information for cases D.1 to D.5.The conventional PID method tracks the desired tip angle with a more significant margin of error than the FPID algorithm.Moreover, it is observed that tracking relatively higher frequency trajectories is challenging for both FPID and PID methods.In contrast, the FFIT2FOF-PID algorithm tracks the higher frequency trajectories accurately, which is primary attributed to its IT2FL membership functions and feedforward compensator.
Moreover, the control efforts of the FFIT2FOF-PID method are notably smaller than those of the FPID and PID methods for all trajectories considered.Both the FPID and PID controllers exhibit noticeable chattering in the control voltage, while the FFIT2FOF-PID method yields smooth variation in the control effort for all considered cases.This is mainly attributed to the fractional operators integrated within the FFIT2FIT-PID controller that help significantly alleviate the chattering caused by the noise in the image sensor signal.
The relative tracking performance of the FPID, PID, and FFIT2FIT-PID control methods are further quantified in table 6 in terms of IAE, IATE, and CE indicators.All of the   performance indicators were evaluated over the entire 30 s test duration.The results also confirm the performance superiority of FFIT2FIT-PID method over the FPID and PID methods for all trajectory cases considered.For instance, for the D.1 case, the IAE indicator for the FFIT2FIT-PID method is 15.2% and 28.7% lower than those of the FPID and PID methods, respectively.The FFIT2FIT-PID method, in-particular, yields considerably better performance for higher frequency trajectory variations.For instance, in the case of trajectory variation at 2 Hz (case D.3), the ITAE index of the FFIT2FIT-PID shows tracking enhancements of about 78.8% and 87.9% when compared to the FPID and PID control methods, respectively.The results summarized in table 6 also show similar trends in the CE index for all the cases considered.

Conclusion
This paper aimed at developing a time-dependent model with a low computational effort and high simulation speed capable of capturing the nonlinear behavior of an MSCR, in addition to designing a robust control strategy that is favorable for tracking relatively high-frequencies trajectories with significant accuracy.In this regard, an accurate FE model was initially developed to estimate the nonlinear behavior of an MSCR under a uniform magnetic field perpendicular to the robot's length ranging from 0.5 − 32 mT.The comparison study between simulation results and experimental findings revealed a maximum deviation of 2.7%.An equivalent model, which benefited from a combination of weighted reducedorder linear models, was subsequently formulated to estimate the robot's nonlinear behavior efficiently.The developed model at the worst-case scenario consisting of 36 linear models with the order of six, which estimated the MSCR's behavior with high accuracy, was implemented in about 6.4 s which was 792 times faster than the FE simulation with a similar simulation time.This favorable advantage was used to design the FFIT2FOF-PID control algorithm, whose control parameters are tuned by GMCOA.The HiL experimental study on predefined trajectories revealed the superiority of the designed control algorithm over the FPID and PID methods, particularly at relatively high-frequencies where those control methods exhibited poor control performance.The developed model and control methodology can provide essential guidance for practical navigation of MCSR in a confined and unstructured environment with uncertainty.The development of a modelbased control strategy based on the developed time-dependent model in environments with uncertainties and disturbances may be considered as a future research study.

Figure 1 .
Figure 1.(a) Minimally invasive treatments, (b) schematic of the distal portion of an MSCR consisting of a soft texture polymer and NdFeB hard magnetic particles, and (c) navigation of an MSCR in blood vessels illustrating bending of the distal portion under an applied magnetic field.

Figure 2 .
Figure 2. A trained trajectory in the 2D state-space with different operating points.

Figure 3 .
Figure 3.The schematic diagram of the interval type-2 fuzzy logic network.

Figure 5 .
Figure 5. Illustration of (a) the mold designed in SOLIDWORKS software, (b) the mold fabricated using an SLA 3D printer, and (c) the fabricated MSCR sample.

Figure 7 .
Figure 7. (a) Scanning electron microscope image of the sample with 17% volume fraction of NdFeB particles, and (b) mechanical tensile test results and Neo-Hookean fitting model.

Figure 8 .
Figure 8. Hardware-in-the-loop experiment design for real-time actuation of the MSCR.

Figure 9 .
Figure 9. (a) FE model of the MSCR and its surrounding airbox, and (b) The transverse angular deformation of the MSCR's tip under different magnitudes of uniform magnetic flux density (peak angular deformation is also shown in the legend).

Figure 10 .
Figure 10.The images of the steady-state transverse angle response of the MSCR's tip (uniform magnetic flux density = 2−25 mT).

Figure 11 .
Figure 11.Comparison of the peak tip deflection response predicted from the model with the measured data under different magnitudes of magnetic flux density.

Figure 13 .
Figure 13.(a) The best cost value among the five times running the GMCOA algorithm, and (b) IT2FL surface rule after optimization with GMCOA.

Figure 14 .
Figure 14.IT2FL network characteristics after optimization with GMCOA: (a) MFs for first input, (b) MFs for second input.

Table 2 .
Geometrical and mechanical specifications of the MSCR sample.

Table 3 .
Comparison of transverse tip deflection response of the FE model with the measured data under different magnitudes of magnetic flux density.

Table 4 .
Reported RSME and STD indicators for developed FE and TPWL models.

Table 6 .
Relative performance indicators for the PID, FPID and FFIT2FOF-PID closed-loop system of the MSCR, evaluated from the experiments.