Task Accuracy Enhancement for a Surgical Macro-Micro Manipulator With Probabilistic Neural Networks and Uncertainty Minimization

Accurate robot kinematic modelling is a major component for autonomous robot control to guarantee safety and precision during task execution. In surgical robotics complex robotic structures and actuation mechanisms are generally employed, therefore machine learning techniques can be adopted to build the model of the robot. Probabilistic neural networks are a class of learning approaches that provide information about the uncertainty of the learnt models. In this work we compare two different probabilistic neural networks (Bayesian and Evidential Neural Networks) to model the kinematics of a surgical robotic instrument and propose a control strategy based on Hierarchical Quadratic Programming (HQP) capable of exploiting the model uncertainty to improve the accuracy and safety of the controller. Simulation and real world experiments on different autonomous path tracking tasks show that the model uncertainty highly affects the control performances and prove the effectiveness of the proposed controller in improving task execution. Note to Practitioners—The push towards reducing invasiveness and patient’s traumas in surgery has lead to the requirement of miniaturized and highly articulated robots. This however comes at the cost of having systems that are hard to model and control, which is one of the major limitations for autonomy in surgical robotics. Machine learning has become very effective in modelling complex systems and probabilistic approaches additionally allow estimating the confidence of the learnt model. In robotics field where high precision is required to perform an autonomous task, like in minimally invasive surgery, the robot model needs to be very accurate and controllers need to guarantee safety in performing the desired task, while satisfying additional motion constraints imposed by the application scenario. This work proposes the use of probabilistic neural networks to model the complexity of a surgical robotic instrument and a control strategy capable of ensuring safety by maximizing model’s confidence and guaranteeing satisfaction of imposed motion constraints. In this work a macro-micro manipulator setup is employed, consisting of an articulated surgical robotic instrument connected to a serial-link manipulator. The proposed modelling and control approaches can be used in any other field where controllers need to highly rely on the robot model due to limitations in using external sensors and where leveraging information about model’s confidence can be beneficial. Currently, the work focuses only on pure kinematic modelling and control, thus neglecting any possible interaction with the environment. Future work will focus on addressing this limitation in order to ensure proper force control and effective autonomy.

Abstract-Accurate robot kinematic modelling is a major component for autonomous robot control to guarantee safety and precision during task execution.In surgical robotics complex robotic structures and actuation mechanisms are generally employed, therefore machine learning techniques can be adopted to build the model of the robot.Probabilistic neural networks are a class of learning approaches that provide information about the uncertainty of the learnt models.In this work we compare two different probabilistic neural networks (Bayesian and Evidential Neural Networks) to model the kinematics of a surgical robotic instrument and propose a control strategy based on Hierarchical Quadratic Programming (HQP) capable of exploiting the model uncertainty to improve the accuracy and safety of the controller.Simulation and real world experiments on different autonomous path tracking tasks show that the model uncertainty highly affects the control performances and prove the effectiveness of the proposed controller in improving task execution.
Note to Practitioners-The push towards reducing invasiveness and patient's traumas in surgery has lead to the requirement of miniaturized and highly articulated robots.This however comes at the cost of having systems that are hard to model and control, which is one of the major limitations for autonomy in surgical robotics.Machine learning has become very effective in modelling complex systems and probabilistic approaches additionally allow estimating the confidence of the learnt model.In robotics field where high precision is required to perform an autonomous task, like in minimally invasive surgery, the robot model needs to be very accurate and controllers need to guarantee safety in performing the desired task, while satisfying additional motion constraints imposed by the application scenario.This work proposes the use of probabilistic neural networks to model the complexity of a surgical robotic instrument and a control strategy capable of ensuring safety by maximizing model's confidence and guaranteeing satisfaction of imposed motion constraints.In this work a macro-micro manipulator setup is employed, consisting of an articulated surgical robotic instrument connected to a serial-link manipulator.The proposed modelling and control approaches can be used in any other field where controllers need to highly rely on the robot model due to limitations in using external sensors and where leveraging information about model's confidence can be beneficial.Currently, the work focuses only on pure kinematic modelling and control, thus neglecting any possible interaction with the environment.Future work will focus on addressing this limitation in order to ensure proper force control and effective autonomy.
Index Terms-Probabilistic neural networks, redundant robots, optimal control, model uncertainty, minimally invasive robotic surgery.

M INIMALLY invasive surgery (MIS) has been replacing
open surgery in many clinical procedures, since it has been shown to cause less bleeding, reduce size of incisions and patient's trauma, and shorten recovery time [1], [2].However, surgical tools are inserted into the patient's body through small incisions, leading to loss of vision for the surgeon and less intuitiveness in manipulating the instruments.
Robotic Assisted Minimally Invasive Surgery (RAMIS) was then developed to overcome these issues and make it easier for surgeons to perform surgical operations thanks to robots' ability to achieve higher precision, higher dexterity, and remove human tremors.Furthermore, employing robots in a clinical environment allows for increased autonomy, which, in turns, has the benefits of simplifying surgical procedures, reducing the workloads for surgeons to avoid fatigue, as well as improving task accuracy and precision [3].Nevertheless, one of the major limitations on the widespread use of autonomous surgical robots is the lack of precise and robust modelling and control techniques needed to ensure safety and precision.
In recent years, a major translation has been towards flexible surgical robots, in order to further overcome the limitations of MIS and augment surgeons' capabilities.The advancements in design and manufacturing have enabled the possibility to build highly complex structures such as highly articulated robots [4], continuum robots [5], and soft robots [6].These structures use complex actuation mechanisms such as tendon sheath, wire drives, or customized joint designs which make them very dexterous, flexible, and well-fitted for navigation in confined areas of the human body during MIS.Despite the promising capabilities of these systems, their applications in real clinical scenarios are still limited due to many limitations, such as lack of effective motion control strategies and proper modelling caused by the high complexity of their structures [7], [8].
From the control side, the control strategies should be able to guarantee accurate motion execution while satisfying various constraints imposed by the surgical procedures [9].In RAMIS, the small incisions where the surgical instrument is inserted restrict the motion of the robot, which is not allowed to move tangentially to the hole plane.Yet, it can only pivot and translate in the insertion direction about the insertion point.One of the challenges for the control of surgical robots is guaranteeing safe execution of the constrained motion at the insertion point, also known as Remote Center of Motion (RCM) [10].This constraint leads the robotic system to lose two Degrees of Freedom (DOFs) [11].A common approach is to use traditional serial-link manipulators on which the surgical instrument is mounted [12], [13], [14], creating a macro-micro manipulator.By exploiting the kinematic model of the macro-micro manipulator, appropriate control strategies can be designed to fulfill the desired task.
However, the complex structures employed in surgical robots make system modelling very challenging.Machine learning techniques have become very popular and effective in dealing with complex physical models, and they have been widely used in robotics for modelling and control [15].
Model-free approaches, such as reinforcement learning, do not require any knowledge about the robot model but instead learn control policies based on some reward function [16], [17].Compared to supervised learning, these approaches are less data-efficient, they require an appropriate reward function definition, and may need to be retrained for new tasks [18].
On the other hand, model-based techniques are used to generate an approximation of the robot model, without the need of analytical mapping, which may be hard to obtain due to the complexity of the system.The model can then be employed to formulate the control as an optimization problem [19] and exploit the efficacy of standard and consolidated control techniques.
What is more, an accurate robot model is highly needed in fields like surgery in order to ensure proper surgical tool tracking and scene reconstruction.Generally, endoscopic cameras are present and they have been widely used for scene segmentation and for tracking the surgical tools [20], [21].However, these estimations do not work in case of reflections due to the tool material, smokes in the environment (that might be caused by the surgical procedure), in case of occlusions due to tissues or organs, or even self-occlusions by the robotic systems themselves.Researchers have been trying to solve these issues by integrating camera-based approaches to kinematic models [22], [23].Having a good kinematic model is therefore necessary to enhance scene reconstruction approaches.

A. Related Work
Different works have focused on model learning for robotics, with the vast majority focusing on learning the inverse kinematics [24], [25], [26], directly computing the control variables from task space variables.Learning the inverse kinematics is a challenging task as there could exist multiple solutions [15], it does not allow the exploitation of redundancy in highly articulated robots [27], and it requires information about the robot pose, which might be unavailable during use due to the lack of sensors, as in the case of robotic surgery.
Learning the forward kinematics of the robot, instead, can be beneficial to overcome these limitations.Additionally, a forward kinematic model allows integration with other robot models, such as in a macro-micro manipulator setup [28].In [30] and [31] Feedforward Artificial Neural Networks (ANNs) were used to learn the forward kinematics of a robotic system and analytical derivation is carried out to compute the robot Jacobian for control purposes.A similar approach is used in [32] and [33] where the derivation of the Jacobian is exploited to solve the redundancy problem by means of Null Space Projection [33].In [35] and [36], instead, the Jacobian computation is directly included during model training.
The main drawback of these machine learning approaches is that they are deterministic, and, as such, they provide only a punctual estimation of the model.However, two kind of uncertainties exist in model learning that refer to noise in the data (aleatoric uncertainties) and uncertainties in the model itself (epistemic uncertainties).Probabilistic models can be used to express these uncertainties, which can then be incorporated into planning and policy evaluation [36].One of the most popular technique is Gaussian Process Regression (GPR) [37], which is a nonparametric regression approach that relies on assuming a Gaussian distribution over functions.The output of the GPR is a Gaussian probability distribution whose variance provides an estimation of the model's uncertainty.
The major drawback of GPR is being very computationally expensive, as it requires inverting the covariance matrix, whose size equals the number of data samples used for training.Bayesian Neural Networks (BNN), instead, are similar to ANN but use a Bayesian regression approach for learning by assigning a prior probability distribution over the network's weights [38].They thus allow building complex mappings but are less prone to overfitting and, most importantly, return information about the model uncertainty.
Recently, Amini et al. [39] have proposed a novel deep learning probabilistic approach for regression, named Evidential Neural Networks (EvNN), that does not make any assumption either on the data distribution, nor on the network's weights.Instead of placing priors on network's weights, EvNN place priors directly over the data likelihood function.By training a neural network to output the hyperparameters of the higher-order evidential distribution, a grounded representation of both epistemic and aleatoric uncertainties can then be learned without the need for sampling, which is a major limitation for BNN.
Being capable of learning accurate models and exploiting the model unceratinties has been shown to be beneficial for control purposes [40], [41] in order to achieve higher accuracy and stability.This is of uttermost importance in fields like RAMIS, where as little damage to the patient as possible should be caused.

B. Contributions
The purpose of this manuscript is to propose a control strategy for a macro-micro manipulator, consisting of a serial-link manipulator and a surgical robotic instrument, capable of exploiting model uncertainties in order to improve autonomous task execution.Since the kinematics of serial-link manipulators can be considered accurate, in this work only the kinematic model of the Micro-IGES [42] surgical robot is learnt by means of BNN and EvNN.The two approaches have been chosen over GPR due to their greater novelty and being less computationally expensive, which is needed to ensure realtime control.The contributions of this work are therefore: • employing and comparing BNN and EvNN for robot kinematic modelling; • integrating the learnt model of a surgical robotic tool with a KUKA serial-link manipulator in macro-micro manipulator structure; • proposing a control strategy based on Hierarchical Quadratic Programming (HQP) to control hyperredundant macro-micro manipulators to perform a desired motion task in full autonomy, while subject to the RCM constraint, and exploiting the model uncertainty to achieve better performance.In [40] a similar approach based on Hierarchical Model Predictive Control (Hi-MPC) with BNN to model and control the Micro-IGES surgical robot was proposed.In this current work we are advancing our research in different ways.Firstly, the Micro-IGES model is here learnt by exploring as much as possible its whole workspace, whereas in [40] the model was learnt on task-specific data and the controller tested on similar paths to those used for model learning.Secondly, here the learnt model of the surgical robot is used in a macro-micro manipulator setup.The Micro-IGES model is however learnt independently of the KUKA, since the kinematic model of the KUKA can be easily obtained from standard approaches.To the best of the authors' knowledge, there is no work in the literature focusing on controlling a macro-micro manipulator combining the kinematic model of a serial-link arm and the learnt model of the attached surgical robot.From a control perspective, this combination results in a larger number of control parameters and larger number of prioritized tasks to include.The employed Hi-MPC would therefore need an additional hierarchical layer, thus increasing its computational cost.HQP is instead faster since it allows solving the control problem directly at the velocity level by means of the Cartesian Jacobians and does not consider the evolution of the system dynamics over time, even though the optimality of the solution might be reduced.Finally, in this work we employ both EvNN and BNN for model learning.EvNN for regression have only been recently presented and have never been used in robot modelling and control.Besides, a comparison between EvNN and BNN for robot modelling is still not present in the literature.
Comparing to the recent work [43] where HQP control has been used on a surgical macro-micro manipulator, this current work adds an additional layer in the subtasks for uncertainty minimization.Furthermore, in [43] an optimization framework for finding the optimal counting connection between the serial-link arm and the surgical instrument is proposed, but the surgical robot's model is not learned.
The manuscript is thus structured as follows.Section II presents and compares the two probabilistic approaches BNN and EvNN for regression, along with a summary of the most widely used kinematic control strategies for redundancy resolution in robotics.Section III describes the kinematic model of the macro-micro manipulator.Section IV introduces the proposed control strategy with uncertainty minimization.Experimental results, carried out both in simulation and on the real robot, are presented in Section V. Finally, conclusions are drawn in Section VI.

II. ROBOT KINEMATIC CONTROL AND PROBABILISTIC NEURAL NETWORKS
In this section, the most common control strategies used in robot kinematic control are presented, along with Bayesian and Evidential Neural Networks for regression.

A. Redundancy Resolution
Given a desired Cartesian motion task h ∈ R m and the robot's forward kinematic model F K (q) mapping the robot's joints q ∈ R n j to the Cartesian space, the goal of inverse kinematics is to compute the proper joint angles q, velocities q, and, eventually, accelerations q that allow the system to execute the desired motion task.The linearity at the velocity level is usually exploited to compute joint velocities from: where J is the Jacobian matrix and ḣ ∈ R m is a desired velocity task.
In case of redundant robots, the number of available joints n j to execute the trajectory is larger than the number of constraints imposed by the desired task m.Therefore, the solution to problem (1), if it exits, is not unique.Different techniques have been developed to solve this problem, such as the Extended Jacobian technique [44] or the Augmented Jacobian [45], but they all place different subtasks at the same level, without any consideration for conflicting subtasks.
If prioritization of subtasks is desired, such that subtasks with higher priority are solved first and then the others are solved with the exceeding DOFs, one common approach consists in using the Gradient Projection Method (GPM) [46].
Recently, especially in the field of humanoid robotics, where researchers have to deal with high level of redundancy, a common approach to solve stacks of prioritized subtasks subject to different constraints (both equality and inequality) consists in applying the Hierarchical Quadratic Programming (HQP) framework introduced by [47], where the control problem is formulated as: where N is the number of prioritized subtasks, λ is a constant coefficient for regularization, ḣn is a desired subtask to achieve, A in,n , b in,n describe additional inequality and equality constraints for each motion subtask.A in,0 , b in,0 , J h,0 can be set to zero or empty, as only used for initialization.This formulation allows finding the optimal solution q * n , without violating the higher priority subtasks 0 . . .n − 1 and it is a generalization of GPM in the presence of equality and inequality constraints.

B. Bayesian Neural Networks
Similarly to Feedforward Artificial Neural Networks (ANN), BNN are capable of representing complicated behaviours, without the need of knowing any mathematical or physical model.ANN, however, are highly influenced by outliers in the dataset [48], suffer from overfitting, and it is difficult to control their complexity [49].BNN, instead, allow including uncertainties in the model, by adding priors to the weights of the network.BNN use probability distributions to quantify uncertainty in inferences and output a probability distribution expressing the belief on how likely the different predictions are [49].
Given a dataset of D observations D = { ỹ1 , . . ., ỹD } and the corresponding input vectors x 1 , . . ., x D , the resulting neural network model is defined as ŷ = y(x, w), where w is the vector of network's weights and biases.In the Bayesian learning, a Gaussian conditional distribution is assigned to each data point p( ỹ|x, w, β) = N ( ỹ|y(x, w), β −1 ), with the output of the network as mean and inverse variance β.Assuming each point independently, the overall distribution on the dataset is computed as p(D|w, β) = D d=1 N ( ỹd |y(x d , w), β −1 ).A prior distribution p(w|α) = N (w|0, α −1 I ), with zero mean and inverse variance α, is also assigned to the network's weights.The posterior distribution on the weights is then: which is non Gaussian, due to the nonlinear network model y(x, w).However, a Gaussian approximation can be found by using different techniques, such as Laplace approximation or Variational Inference [50].In this work, we adopted the latter one, which allows identifying the Gaussian approximation of (3) as: q(w|D) = N (w|w * , A −1 ), with the optimal weights w * as mean and A as inverse covariance matrix, utilizing Kullback-Leibler divergence.This is performed during offline training by employing the Dense Flipout estimator [51] for each network layer.
Once the Gaussian approximation for the weights' posterior is computed, the predictive marginal likelihood on the data can be obtained as p( ỹ|x, D) = p( ỹ|x, w, β)q(w|D)dw.However, the analytical solution is intractable due to the nonlinearities in the network mapping.One way to solve this issue is by using Monte Carlo sampling, which, however, might be very computationally expensive as multiple samples are needed.Instead, by linearizing the network mapping around the optimal weights w * , the predictive distribution on the network output can be computed analytically as [50]: where g w = ∂ y(x,w) ∂w is the gradient of the network over the weights.The output of the BNN is therefore a Gaussian distribution with the network model ŷ = y(x, w * ) as mean, and σ 2 N N as variance.The variance results to consist of two terms: the aleatoric uncertainty of the data σ 2 al = β −1 and the model epistemic uncertainty σ 2 ep = g T w A −1 g w .Thanks to the analytical expressions, it is possible to compute the network derivatives in a closed form, both for the prediction and the uncertainty σ 2 N N , as shown in Appendix A.

C. Evidential Neural Networks
Even though BNN allow computing an analytical closed-form expression for the predictive distribution, it relies on different assumptions on the prior data and weights' distribution, and on some approximations to make the analytical derivations tractable.In order to overcome these issues, recently Amini et al. [39] have proposed a new probabilistic neural network framework named Evidential Neural Networks (EvNN).
EvNN assume that the data follows a Gaussian distribution, but with unknown mean and variance, i.e. p( ỹ|x, w, μ, σ ) = N ( ỹ|μ, σ 2 N N ).The mean μ and variance σ 2 N N are assumed to be drawn from a Gaussian distribution Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I COMPARATIVE SUMMARY OF BNN AND EVNN
and an Inverse Gamma distribution respectively: where (•) is the gamma function, α > 1, β,υ > 0 and γ ∈ R are parameters which are learnt along with the network's weights.
The posterior distribution of the parameters μ, σ 2 N N is then approximated as: ) , (6) which is the Normal Inverse-Gamma (NIG) distribution.Drawing a sample μ j , σ 2 N N, j from the NIG distribution returns one instance of the data likelihood N ( ỹ|μ, σ 2 N N ) and the hyperparameters γ , α, β, υ determine the location and the uncertainty associated with the inferred likelihood.
Similarly to BNN, EvNN consist of a cascade of feedforward layers, with the addition of a final layer called Evidential Layer.This layer maps the output of the feedforward networks onto a 4 dimensional vector returning the quantities γ , α, β, υ.
Finally, given the NIG distribution it is possible to easily estimate the model prediction ŷ, aleatoric uncertainty σ 2 al , and the model epistemic uncertainty σ 2 ep as: Being all the quantities computed directly from a feedforward layer structure, an analytical expression can be obtained, with consequent easy computation of the derivatives (Appendix A).

D. Network Comparison
Table I reports a summarizing comparison of the two probabilistic approaches, BNN and EvNN.The main difference between the two approaches is that BNN assume a known distribution for the data likelihood and weights prior.In EvNN, instead, the data is assumed to be Gaussian, but its parameters are unknown.
Both frameworks are built on a feedforward network layer structure and this allows for a closed-form analytical expressions for the predictive distribution mean and variance.However, the computational cost (especially for the derivatives) in BNN is larger than in EvNN.
With regards to network training, BNN require complex techniques such as variational inference, Flipout, or Monte Carlo sampling, due to the intractability of computing the weights posterior analytically.In EvNN, instead, the network is simply trained using a traditional regression cost function defined in [39].
Finally, in terms of output prediction, BNN lead to generally more conservative results, with smaller changes in the epistemic uncertainty.EvNN, on the other hand, have larger variations and, thus, manage to explain better the model uncertainties.

III. ROBOTIC SYSTEM DESCRIPTION
In this section an overview of the robotic system setup is presented, describing the Micro-IGES surgical robotic tool and the KUKA LBR IIWA robotic arm.

A. Micro-IGES Surgical Robotic Tool 1) Robot Description:
The Micro-IGES is a custom-made surgical robotic tool, composed of a rigid shaft and a flexible section.In total it has 6 DOF for motion control and 1 gripping DOF, with the two elbow joints consisting of a pair of coupling joints.The coupling of the two pairs of joints of the elbows occurs at the distal driving unit at the end of the shaft.As shown in Figure 2, the robot's base is attached to the motor package by means of pinned connectors on the capstans, and the tendons are routed on the capstans themselves.Due to the current setup, in this work only 5 DOF are considered (Roll, Elbow 1, Elbow 2, Wrist Pitch, Wrist Yaw).
2) Robot Kinematic Model: Because of the generally nonlinear motor to joint (and joint to motor) mapping q u = f(θ u ), being θ u ∈ R n m the vector of motor positions, P θ u ∈ R n m the vector of motor velocities, and q u ∈ R n j the vector of joint positions, the kinematic model of the robot, with respect to its base frame {R F b }, can be rewritten as: where b T q , b T u ∈ R 4×4 are the Cartesian end-effector pose as a function of the joint and motor values respectively, b v u ∈ R 6 the Cartesian end-effector twist (linear and angular velocities vector), b J q is the Cartesian task Jacobian with respect to the joint variables, and b J u the Jacobian with respect to the controllable motor values.The matrix L, with L i, j = ∂q u,i ∂θ u, j , is the motor to joint differential matrix [52].The control problem for the Micro-IGES can therefore be formulated as a function of the motor values, which can be directly measured and controlled.As described in III-A1 the system state can is expressed by Roll, Elbow 1, Elbow 2, Wrist Pitch, Wrist Yaw.
3) Learning-Based Kinematic Modelling: Due to the high nonlinearities in the system caused by the tendon-transmission, an analytical model of the robot is hard to compute.In [53] tendon nonlinearities are compensated using friction models like LuGre or Wouc-Ben, yet, the high number of joints and the lack of external sensors for the online estimation of the joint configuration for the Micro-IGES robot make such an approach hard to apply.Even though a simple backlash compensation strategy is implemented in the robot (whose description is out of the scope of this work), a proper mapping from actuation, to joint, to Cartesian space is still required.Since the tip position can be easily measured, BNN and EvNN are employed in this work to build the mapping from motor values θ u to 3D Cartesian tip position b Pu,b ∈ R 3 , with respect to the Micro-IGES base frame.By modelling the robot kinematic model offline, external sensors are not needed of online operations.
For both BNN and EvNN three independent networks with the same input (the Micro-IGES motor values) have been used to learn the kinematics on each Cartesian component independently.The output of each network is therefore a one-dimensional vector corresponding to each Cartesian component.Building the forward kinematic mapping by means of BNN and EvNN also allows estimating the robot Jacobian . Each network provides one single row of the estimated Jacobian.Additionally, BNN and EvNN models output the uncertainty on the estimated tip position σ (θ u ) ∈ R 3 .Thanks to its analytical expression, it is possible to compute its derivative with respect to the control variables J σ,u = ∂σ 2 ∂θ u ∈ R 3×5 .In this work we are focusing only on learning the mapping to the Micro-IGES's tip position for two main reasons.This is because learning the orientation needs further investigation, especially because the network architecture needs to be designed in such a way that the derivatives of the output are consistent with kinematic quantities used for control purposes like the end-effector's angular velocity.A parallel research is being conducted to address this limitation.

B. KUKA Robotic Arm
The KUKA LBR IIWA is an articulated industrial robot with 7 joints.Thanks to its simple design and direct transmission (the motors are placed directly at the joints), the kinematic model of this robot can be easily computed by means of standard kinematic approaches (e.g.Denavit-Hartenberg convention).Therefore, the kinematics of the KUKA can be expressed with respect to its base frame {R F 0 } as: where q k ∈ R 7 is the KUKA joint position vector.

C. Macro-Micro Manipulator Kinematics
The macro-micro manipulator is composed by the Micro-IGES and the KUKA robots (Figure 1, 3), with the Micro-IGES motor package directly attached to the KUKA's end-effector.In total, the system has 12 DOF (7 for the KUKA and 5 for the Micro-IGES) and its state is then described by q = q k θ u T .
By expressing the kinematics models of the serial-link arm and surgical robot with respect to a common reference frame (e.g. the KUKA base frame), it is possible to obtain the macromicro manipulator's end-effector position as [43]: where 0 P k is the KUKA's end-effector position in its base frame, k P b,k is the Micro-IGES' base position with respect to the KUKA's end-effector frame, b Pu,b is the estimated Micro-IGES' end-effector position with respect to its base frame resulting from the neural network models, 0 R k is the KUKA's end-effector rotation matrix in its base frame, and k R b the Micro-IGES' base rotation matrix in KUKA's endeffector frame.
By deriving (10) with respect to time and with some mathematical manipulations [43], it results that: where J ω,i for i = 1 . . .7 are the KUKA's end-effector angular Jacobians for each link, J v,k and Ĵu the linear Jacobians of the KUKA and the Micro-IGES respectively.Equation ( 11) allows for simple computation of the Cartesian Jacobian of the combined KUKA and Micro-IGES Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Ĵ u is the estimated Micro-IGES Jacobian from the neural network, expressed in the KUKA base frame.
Ĵ K U K A ∈ R 3×7 and Ĵ u Iges ∈ R 3×5 describe the contribution of each robot to the end-effector's motion.

IV. CONTROL METHOD
In this section, the control strategy to execute a desired task by using the combined KUKA and Micro-IGES robots is presented.

A. Stack of Tasks
In order to optimally exploit the hyper-redundancy of the macro-micro manipulator, different motion subtasks can be specified.Similarly to [43], here we employ the HQP controller (2) since a strict hierarchical prioritization ensures that the primary subtask is always executed, while avoiding the need of manually updating the weights whenever a new motion is specified (which instead affects soft task prioritization schemes [54]).
In a surgical application like tumor resection at least two motion subtasks can be specified: keeping the Remote Center of Motion (RCM) fixed and following a desired path (in order to remove the tumor).

1) Remote Center of Motion:
The RCM is a point, usually the surgical tool insertion point, where only rotations and insertion motion can be achieved.During laparoscopic procedures, for instance, the RCM coincides with the location of the hole on the patient's body where the surgical tool is then inserted as in Figure 4. Therefore, the motion at the RCM needs to be constrained along the surface tangential directions (y, z in Figure 4).
The Cartesian position and velocity of the RCM with respect to the robot base is computed similarly as position with respect to the KUKA's end-effector, in the end-effector's frame and 0 J RC M a matirx expressing the contribution of the KUKA to the RCM.
In [55] the RCM motion is expressed with respect to the serial-link manipulator's end-effector frame.However, due to the motion of the KUKA arm, in this work the RCM motion is expressed with respect to the insertion point frame in Figure 4. Therefore the rotation matrix of the hole frame with respect to the KUKA base frame.Thanks to this transformation, assuming the hole fixed, the RCM motion to consider is only along the y, z components.Consequently, h ḣRCM = 0 0 T in order to guarantee the RCM motion in the tangential directions to be fixed and h J RC M ∈ R 2×12 .In order to constrain the RCM, the primary subtask cost function can be expressed as c RC M = || h ḣRCM − h J RC M q|| 2 .However, in cases where the RCM is not stationary (like patient breathing) it is still possible to track the RCM by properly specifying the velocity h ḣRCM .In these cases the RCM velocity needs to be known either a priori or by measurements with external sensors.
2) Path Tracking: The secondary motion subtask is executing a desired motion, which, in cases like tumor resection, could be having the macro-micro manipulator's endeffector follow a desired path.Given reference Cartesian path hp ∈ R 3 , and its the corresponding desired velocity ḣp , from (11), the cost function associated to the secondary motion task can be expressed as c p = || ḣp − Ĵ q|| 2 .

B. Prioritized Motion Control
In order to guarantee accurate execution of the two prioritized subtasks the HQP problem ( 2) is modified as follows: where n = 1 corresponds to guaranteeing RCM motion ( ḣ1 = h ḣRCM and J h,1 = h J RC M ), and n = 2 is for the path tracking ( ḣ2 = ḣP and J h,2 = Ĵ).Even though the RCM task could be added as a constraint for the secondary task, the HQP formulation favours solving the prioritized tasks in a least square sense in order to avoid conflicts with the constraints of the subsequent subtasks [56].For the first subtask, J 0 can be set to zero or empty as only used for initialization.q m , q M are the joints bounds, and dt the motion sampling time.The matrices W t,n , W K U K A,n , W u Iges,n are diagonal weighting matrices which are used for different purposes.
For each subtask n, they can be assigned different values with W t,n ∈ R 3×3 assigning different importance to each task dimension, whereas W K U K A,n , W u Iges,n ∈ R 3×3 allow penalizing the contribution to the end-effector's motion of the whole robotic structure and they can be used to make only one single robot move.Currently, the entries of the weighting matrices are manually set based on an heuristic trial and error approach.Optimal ways to assign the weights are however being explored to better automate and improve the process.Additional constraints, for instance, on the joint velocities could also be included, yet an appropriate reformulation would be needed in order to ensure there are no conflicts in the kinematic constraints [57].

C. Uncertainty Minimization
In a surgical environment precisely tracking a trajectory is fundamental to guarantee safety of the patient and reduce traumas.Yet, the tracking accuracy depends on the precision and reliability of the kinematic model.Thanks to the use of the two probabilistic neural network models (BNN and EvNN), it is possible to estimate the uncertainties in the kinematic model and control the robot to move within safe state-space regions [58].For that purpose, the additional task of minimizing the kinematic model epistemic uncertainties σ can be added and solved as a lower priority subtask for the path tracking.This means that the uncertainty minimization subtask will try to exploit the robot's redundancy and self-motions to reduce the model's uncertainty, without deteriorating the RCM and path tracking subtasks.
The HQP problem ( 12) is solved at each time step t discretizing the motion to be performed.The uncertainty minimization could be directly included in the HQP framework as a third subtask.However, in this case the uncertainty would be minimized only with a single step at each motion time step t.Due to the high nonlinearity of the function itself, a one-step minimization (OneStep) may not be enough.
Therefore, to enforce larger minimization a Sequential Quadratic Programming (SQP) problem is iteratively solved at each time step t.This corresponds to subsampling each time step in smaller steps at which the commanded motion task would not vary, but during which the robot tries to reconfigure itself to a less uncertain state.The optimal solution q * , q * from the HQP with the two prioritized subtasks (RCM and path tracking) will correspond to the end-effector position P * and to the RCM position P * RC M .The uncertainty minimization should not deteriorate the optimal solution from the HQP, therefore, after the initialization q 0 = q * the SQP is formulated at each iteration i = 1 . . .I as: The damping factor λ is used to ensure small changes in the control variables at each iteration, the constraints (13b)-(13c) guarantee that at each iteration both the tip position and RCM position do not change, thus preventing the the optimal solutions from HQP ( P * , P * RC M ) from being altered and causing unwanted Cartesian motions, and (13d) to make sure the limits on the control variables are not violated.In this work the maximum number of iteration I is heuristically set to 100, as a compromise between properly minimizing the model uncertainty and reducing the computational cost.Since only the Micro-IGES model is learnt, only its control variables θ u will affect the model's uncertainty σ .Therefore, , where J σ,u contains the derivatives of the model uncertainty with respect to the Micro-IGES motor positions.
Once the SQP is solved, the optimal solution found is commanded to the robot at time step t. Figure 5 shows a schematics of the control strategy with the SQP for uncertainty minimization.
Since in SQP the model uncertainty σ 2 , its derivative J σ , and the end-effector and RCM Jacobians Ĵ, h J RC M need to be estimated at each iteration i , the computational time tends to be high, especially when BNN is used for modelling.In order to reduce the time, those values are updated only after some iterations (e.g.every 10 iterations).The computational time is however related to the complexity of the network architecture used.We report the computational times for both the EvNN and BNN models in section V-A2 for the architectures used to model our system.
For both the HQP and the SQP problems qpOASES [59] has been employed as solver given its reliability and fast computational speed.

V. RESULTS
In this section, results for the proposed method applied both to simulation and real world tests are presented, comparing the performances of the control strategy when employing both BNN and EvNN for modelling, without and with the proposed uncertainty minimization strategy.

A. Simulation Results
For the simulation tests V-REP simulator [60] has been used and Figure 3 shows the simulation environment.For visualization purposes, the Micro-IGES is directly connected to the KUKA's end-effector, just neglecting a vertical translation.
In the real case, an approximated model of the Micro-IGES can be built by using the same motor to joint mapping as in [61] and by Denavit-Hartenberg (DH) convention to map the joints to the tip position.Due to the many nonlinearities, however, the model is not accurate.In order to make the simulation test as close as possible to the real case, BNN and EvNN are here employed to learn the modelling error between a nominal approximated Micro-IGES model and a ground truth one.The nominal model is generated by using a linear motor to joint mapping; the ground truth model, instead, employs a nonlinear mapping for each motor i of the form q u,i = q u,i,max − q u,i,min where q u,i,max , q u,i,min are the Micro-IGES joint position bounds, and θ u,i,max , θ u,i,min the motor position bounds.Even though a proper mapping from motor to joints is not available for our robotic system, here the function has been chosen in order resemble the results from the analysis performed in [62].
To generate the data, the Micro-IGES motors are excited with sinusoidal motion of the type Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.(sin(2πψ t T ) + 0.2 sin(20πψ t T )) (Figure 6).In order to explore as much as possible of its workspace, each motor is excited with two possible choices ψ = {2, 4}, resulting in a total of 2 5 combinations.An additional independent excitation for each single motor is also included.In total 2 5 + 5 = 37 excitations were commanded and, lasting each motion T = 30 s with a sampling rate dt = 200 ms, 5550 data points were collected.The sampling rate is the same later used for control, and its choice was dictated by the real world experimental setup, as explained in V-B.
While commanding the motor values, the tip positions from the ground truth and the nominal models are collected.Gaussian noise with zero mean and 1 mm standard deviation is also added.It is worth remarking that since the neural networks are used to only model the Micro-IGES, the KUKA robot is not controlled in this part, but acts only as a support.

1) Comparison of Neural Network Modelling:
For training the BNN and EvNN models the collected data is randomly split into a training and test set (80%-20%).In order to learn the kinematic model of the Micro-IGES, both for BNN and EvNN three independent networks are employed to map the motor position values to the Cartesian x, y, z components.For both approaches, and for each Cartesian component, the networks' structure consists of an input layer with 5 neurons (taking as input the Micro-IGES motor position values θ u ), and two hidden layers with 30 and 10 neurons each.For BNN, swish activation function is used, whereas sigmoi d is used for EvNN.This is because sigmoi d activation lead to very small variations in the BNN model uncertainty due to the function saturation; the swish activation, instead, leads to too large variations in EvNN.The choice of the architecture was led by the model performance from previous works (e.g.[40]) and by a trade-off between model accuracy and computational cost.A deeper network with more neurons would achieve higher accuracy, but would lead to higher inference time, especially for computing derivatives needed in the controller.
The Root Mean Square Error (RMSE) between the learnt model and the ground truth simulated model resulted to be 0.863 1.061 0.397 mm for x, y, z components for BNN and 1.002 1.068 0.646 mm for EvNN.Both models manage to accurately reconstruct the ground truth model, with BNN performing slightly better than EvNN.
Figure 7 shows the norm of the model uncertainty σ on the collected dataset with both approaches.As expected, BNN result to be more conservative than EvNN, with a mean norm of the uncertainty σ of 1.82 mm, compared to 2.05 mm of EvNN.
With regards to computational cost, on a Ubuntu 16.04 Intel Core i7-4500U machine, a single query of the BNN to compute the network prediction and its derivatives (both for the tip position and the model uncertainty) takes about 2 ms; a single EvNN query, instead, takes 0.6 ms.We would like to stress that even though the models are trained offline and the control performed in pure open-loop without any sensory feedback (as described in the following), it is important to consider computational costs to ensure the desired task is completed in a desired time frame.Additionally, having a model with low computational costs would allow its usage even in closed-loop control or adaptive modelling approaches.

2) Validation of Uncertainty Minimization:
In this test the proposed SQP (13) for uncertainty minimization is validated while performing a desired motion task with the whole macromicro manipulator.Figure 8 shows the results when using SQP and a one-step minimization of the model uncertainty, both for BNN and EvNN.As shown, if the minimization is directly included in the HQP with the OneStep minimization, the optimizer does not lead to large reductions of the model uncertainty.With the proposed SQP, instead, the model uncertainty can be more successfully decreased: when employing the BNN network model there is no reduction of the mean uncertainty value with the OneStep minimization and a 0.02% reduction with the SQP, with respect to the case without any minimization; when employing the EvNN, the OneStep minimization yields a 0.01% reduction and the SQP a 0.6% reduction with respect to the mean uncertainty value for the case without minimization on this specific task.Understanding by how much the uncertainty should be reduced to lead to improvements in the performance is a hard task and it is still under investigation, yet, as discussed in both the simulation and real world tests, minimizing the uncertainty clearly improves the tracking accuracy [40], [58].
In terms of computational time, adding SQP with the uncertainty minimization increases the amount of time needed to compute the control variables.When employing BNN, solving the HQP (thus querying the network output and derivatives and solving the multitask optimization problem) takes about 16 ms, and the solution of the SQP with I = 100 iterations requires about 90 ms; when using EvNN, instead, the HQP and SQP take 2 ms and 30 ms, respectively.
3) Control Tests: In this test the macro-micro manipulator is required to follow different paths in full autonomy, while keeping the RCM fixed at a specified location.We compare the results when using the nominal simulated model of the Micro-IGES, the learnt model by means of BNN and EvNN, and the learnt models with additional uncertainty minimization in the control.The paths consist of three circles of 5, 10, 30 mm radius, two squares of 20, 60 mm side, and two randomly generated polygons with 3 and 5 vertices.Both the circles and the squares are on a fixed plane at a distance of 0.1 m from the insertion point in the insertion direction.In each test, the Micro-IGES robot is inserted through the hole and it first reaches the center of the path, and subsequently starts following its contour.For each path a desired motion time was specified and it was set to 40 s for the circles, 70 s for the squares, 55 s and 85 s for the random paths with 3 and 5 vertices respectively.
It is worth noticing that the control is carried out completely in open-loop, without any sensory feedback about the manipulator's end-effector position.
The circular and square paths were chosen to reside within the Micro-IGES reachable workspace.The random paths, instead, exceed its workspace, thus requiring more motion contribution from the KUKA arm.For this reason, the penalization matrix in ( 12) is set to b W K U K A,2 = di ag(10 −1 , 10 −1 , 0) for the circular and square paths, and b W K U K A,2 = di ag(10 −2 , 10 −2 , 0) for the random paths, whereas W u Iges,2 = 0.The KUKA penalization matrix is here expressed in the Micro-IGES shaft base frame in order to penalize motions in the tangential directions to the insertion point's surface and make sure that the KUKA properly contributes to the translation motion along the insertion direction.For the primary subtask, instead, W K U K A,1 = 0, being the RCM motion affected only by the KUKA arm.
For the uncertainty minimization, the parameters of the SQP were set to λ = 1 for BNN and λ = 10 3 for EvNN, with a maximum number of iteration I = 100.
Figure 9 shows the tracking results on the different paths and Table II reports the tracking errors and the mean uncertainties over the paths.Evidently, using the inaccurate nominal model of the Micro-IGES results in the largest tracking errors.By learning a more accurate model, instead, the tracking errors are reduced.Furthermore, the proposed uncertainty minimization proves effective in ensuring that the robot moves within more-reliable state-space regions where the model was learnt better, thus yielding further improvements in the path tracking, both when using BNN and EvNN.As it can be noticed from Table II, on every path the HQP+SQP results in smaller mean uncertainties both for BNN and EvNN.Slightly worse performances occur in tracking the largest circle and square paths, them being at the boundary of the maximum reachable workspace of the Micro-IGES.In fact, even the mean model's uncertainties on this paths are higher than on the others.
The main disadvantage of the uncertainty minimization is that it leads to slighter larger RCM tracking errors due to the sequential approximations of the network models.However, in each test the RMSE for RCM tracking is always submillimiter.
Comparing the overall results of the tests performed, the mean total tracking error norm (considering both tip and RCM position) for the BNN resulted to be 3.84 mm and 1.83 mm, when employing the control without and with uncertainty minimization, and 3.36 mm and 2.39 mm for the EvNN.The proposed uncertainty minimization thus led to a −52% tracking error reduction for the BNN and −29% reduction for the EvNN.

B. Real World Results
Since EvNN proved to be satisfactory in the simulated control tracking tests, and being computationally more efficient than BNN, for the real world tests we only employed Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.EvNN to model the Micro-IGES and control the macro-micro manipulator.
1) Robot Model: Similarly to the simulation case, EvNN are used to learn the modelling error between an approximated model of the Micro-IGES and the real one.The approximated model is the same as the one used in simulation, based on Denavit-Harteneberg (DH) convention.In order to collect data, the same strategy as in the simulation was carried out, commanding combinations of sinusoidal motions to the motors, but with a longer period T = 40 s.Electromagnetic (EM) trackers have been used to collect the Micro-IGES tip position data (Figure 1), with a total of 7235 points generated, due to the a sampling rate of 200 ms of the sensors.It is worth noticing that during the data acquisition, the KUKA arm is fixed and acts only as a support.The obtained datapoints therefore describe the motion of the Micro-IGES with respect to its base frame.The same 80% − 20% random splitting as in the simulated model learning was used and Table III reports the RMSE on the training and test datasets between the measured tip position and the network's predicted one.Also in this case the networks' structure for each Cartesian component x, y, z consists of one input layer with 5 neurons 2) Control Tests: For the tracking task, the end-effector is required to autonomously follow a circular path of 10 mm in radius and a square path of 20 mm side on a fixed plane at 65 mm in the insertion direction from the hole, while keeping the RCM fixed at the insertion point.In this scenario, the motion time was set to 2 min for the circle and 3.5 min for the square.This choice was mainly dictated by reducing the stress on the system and avoid possible damage from reaching too high speeds.Notice that the motion tasks are expressed with respect to the macro-micro manipulator's base frame coinciding with the KUKA's base frame.To evaluate the tracking accuracy of the macro-micro manipulator's tip and of the RCM, two EM trackers are used, one at the robot's tip and one on the shaft (Figure 1), whose resolution is of 10 −2 mm.The measurements are only used for assessing the tracking accuracy and do not provide any feedback to the controller which acts in open-loop mode.Since the EM trackers measurements are referred with respect to the EM tracker's receiver frame, the sensory data is then mapped onto the macro-micro manipulator's base frame by knowing the relative poses of the receiver with respect to the desired reference frame.
To better evaluate the control performances, four different tests for each path are run.In each test, however, the same control parameters are used, namely b W K U K A,2 = di ag(1, 1, 0) and λ = 10 3 for the SQP, both for the circular and square paths.Figure 10a shows the paths followed by the robot when employing the initial DH model, the EvNN learnt model without uncertainty minimization in the control, and the EvNN learnt model with the additional SQP in the control.Table IV reports the mean tracking error norm for the tip position, while Figure 10b shows a summarizing error plot over the performed tests.
The variability of the results in each test is due to the repeatability capabilities of the Micro-IGES robot and mainly caused by the nonlinearities in the transmission and the inaccuracy in the initial positioning, which is manually executed.In each test, since they were all run in the same conditions, the mean norm of the model uncertainty σ resulted to be 8.59 mm and 0.85 mm on the circular path, without and with uncertainty minimization respectively, and 6.86 mm and 0.95 mm on the square, thus proving the effectiveness of the SQP in maximizing the model's confidence.From the results it is noticeable that minimizing the uncertainty largely improves the tracking accuracy in both cases.On the circular path, EvNN+SQP resulted in a 33% improvement over using the standard DH approach, and 22% improvement on the square path.With regards to the RCM, instead, the mean RCM error norm resulted to be 0.629 mm, 0.695 mm 0.724 mm for the initial DH model, EvNN, and EvNN with uncertainty minimization respectively on the circular path, and 0.588 mm, 0.603 mm 0.844 mm on the square path.Even though by learning the Micro-IGES model the control performance improve, sources of errors due to possible inaccuracies in the calibration of the matrix from the KUKA end-effector to the Micro-IGES base, possible vibrations and deformations of the rod due to the interactions with the tendons passing through it, and not fully compensated nonlinearities still result in a deterioration of the tracking precision.

VI. DISCUSSION AND CONCLUSION
Achieving high accuracy and precision in robotics is of uttermost importance, especially in fields like minimally invasive surgery, where as little as possible damage to the patient should be caused.Even though machine learning approaches allow building complex robot models and thus improve the controller performance, having an estimation of the model uncertainty can be beneficial.Probabilistic neural networks like BNN and EvNN can, therefore, be preferred over traditional ANN, because they provide a probability distribution as output, with the possibility to extract information about the model uncertainty.
Optimal control techniques like HQP can be employed to utilize the redundancy of the robotic system in order to make the system move within more reliable state-space regions, that were better explored during the network training and that, therefore, lead to a more accurate and less uncertain model.The proposed HQP with additional SQP for uncertainty minimization has proved effective in improving the control performances of a macro-micro manipulator in autonomously following desired paths.However, the tracking accuracy highly depends on how well the neural network has learnt the model.In the simulation tests, both BNN and EvNN satisfactorily learn accurate models and perform well in the path tracking.In the real scenario, instead, the model to learn is more complex, resulting in poorer tracking accuracy.
Even if BNN may learn a better model, its main drawback is the computational time required to query the network output and derivatives.EvNN, instead, directly outputting the parameters of the model's probability distribution, require much shorter computational time.Furthermore, EvNN result in generally higher uncertainty intervals, thus favouring the uncertainty minimization.In BNN the damping factor of the SQP needs to be set to small values in order to properly minimize the model uncertainty.This, however, may result in large jumps of the control variables with consequent inaccuracies in the network linearization, and, thus, failure of the controller.
The proposed approach relies on the solving a quadratic program both in the HQP and in the SQP, thus requiring linearizing the learnt models.In order for the linearization to work well, slow motions need to be commanded, which is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the case in applications like surgery.Nevertheless, nonlinear solvers could be used to achieve better performances, even though the computational cost may increase.In addition, the performances of the controller and of the uncertainty minimization still depend on different parameters, like the robot penalization matrix and the damping factor λ for the SQP.Yet, strategies to optimally choose these parameters can be implemented.To more optimally and effectively minimize the model uncertainty, MPC could be employed since it would consider the evolution of the system in future time steps, and thus perform a similar minimization to the SQP problem but with future predictions.The main issue to address in this case is to have a hierarchical MPC as in [40] but with more layers and a more complex system dynamics, given the inclusion of the KUKA arm.Such an approach is under investigation.Finally, even though our proposed approach shows improvements over using the currently implemented DH model, there are still some nonlinearities in the model that make learning tedious and the control inaccurate.Additionally, in this work both the modelling and control of the whole system have been carried out in free space motion, without any interaction with the environment.Contacts with tissues would most likely change the behaviour of the surgical instrument during its motion.Future work is aiming at tackling this problem by either developing a robot model learning approach that is aware of contacts and that can properly model the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
system behaviour under the identified contacts, or by applying real-time model adaptation techniques.In both cases, however, additional information by external sensors might be needed.

APPENDIX A PROBABILISTIC NETWORK DERIVATIVES
Consider a feedforward neural network with L layers, where y l = a(z l ) ∈ R n o is the output of layer l, a(•) is the activation function, z l = W l x l , with W l ∈ R n o ×n i contains both the weights and biases of the layer, and x l ∈ R n i the layer input.The gradient of the network with respect to all the network's weights w [0:l] up to layer l can be computed as: The derivative ∂ y l ∂ z l ∈ R n o ×n o is a diagonal matrix whose elements depend on the chosen activation function and can be analytically computed.for each element j = 1 . . .n o .
The derivative ∂ z l ∂w l ∈ R n o ×n o n i , for each weight ω l h,k is such that: x l k , if h = j 0, otherwise .
By iterating through all the network layers l = 0 . . .L it is then possible to analytically compute the gradient of the full network with respect to all the network's weights.
The derivatives of the network with respect to the inputs x 0 can be similarly computed analytically and iteratively considering that, for each layer l: Starting again from the first layer l = 0, which is such that ∂ z 0 ∂ x 0 = W 0 , the final network derivative can be obtained.With regards to second partial derivatives, it results that: The partial derivatives ∂ 2 y l ∂ x 0 ∂ z l can be easily obtained analytically for each element j as: where for swish activation, and ∂ z l j ∂ x 0 is the j -th element of W l ∂ z l−1 ∂ x 0 which can be obtained iteratively as in (15).
Given the property of it also turns out that: In addition, considering that 16) can be obtained analytically and iteratively.

A. Bayesian Neural Network Derivatives
In Bayesian Neural Networks the 1D output predictive distribution is computed as: p( ỹ|x, D, α, β) = N ( ỹ| ŷ, σ 2 (x)), and where ŷ is the output of the feedforward network BNN is built upon.The derivative of the network output with respect to the input can then be computed as in (15) as ∂ ŷ ∂ x = g x ∈ R 1×n i0 , with n i0 the dimension of the input of the first layer With regards to the derivatives of the network model uncertainty σ 2 with respect to the input x it results: where g xw ∈ R n w ×n i0 and n w is the total number of weights in the network.

B. Evidential Neural Network Derivatives
EvNN consist of a feedforward network structure, with the addition of an evidential layer as output layer.Given the output of the feedforward network ŷ = ŷ0 ŷ1 ŷ2 ŷ3 T ∈ R 4 , the evidential layer is such that γ = ŷ0 , υ = S P( ŷ1 ), α = S P( ŷ2 ) + 1, and β = S P( ŷ3 ), where S P( * ) = ln (1 + e * ) is the so f t plus function.Due to the higher dimensionality of ŷ, its gradient would be g x ∈ R 4×n i0 .
In EvNN the predictive output is ŷEvN N = γ , therefore its derivative can be computed as ∂ ŷEv N N ∂ x = g x,0 , meaning the first row of the gradient g x .The epistemic uncertainty is instead σ 2 ep = β υ(α−1) which yields: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 1 .
Fig. 1.The macro-micro surgical manipulator, composed by a KUKA arm and the Micro-IGES robot, and the experimental setup.

Fig. 4 .
Fig. 4. Representation of the remote center of motion.

Fig. 5 .
Fig. 5.The control strategy based on HQP for keeping RCM and path tracking, with additional model uncertainty minimization.

Fig. 6 .
Fig.6.Example of the exciting motion for data collection.A is the amplitude of the wave and ψ T is the period. θ

Fig. 7 .
Fig. 7. Norm of the model uncertainty σ on the collected dataset for learnig the simulated Micro-IGES model.

Fig. 8 .
Fig. 8.Comparison of the variance minimization using OneStep optimization and the proposed SQP.

Fig. 10 .
Fig. 10.Path tracking results for the real world tests: (a) the comparison of the paths followed; (b) summarizing box plot of mean tip position error norm || P || in the KUKA base frame.For better visualization, 10a only shows the paths in x, y.
For a sigmoi d function it results that

∂
x the j -th element of g x .

TABLE II RESULTS
FOR THE SIMULATION CONTROL TESTS.|| P ||, || RC M ||, || P || max , ||σ || ARE THE MEAN ERROR NORM BETWEEN THE DESIRED AND ACTUAL ROBOT'S TIP POSITION AND BETWEEN THE DESIRED RCM AND THE ACTUAL RCM POSITION, THE MAXIMUM TIP POSITION ERROR NORM, AND THE MEAN NORM OF THE MODEL UNCERTAINTY.ALL MEASURES ARE IN MM Fig. 9. Path tracking results for the simulation tests, using the initial inaccurate DH model, BNN, and EvNN models.P is the error between the actual and desired robot's tip position and σ the tip position's predicted model's uncertainty.

TABLE III RMSE
IN MM ON THE TRAINING AND TEST DATASETS FOR THE INITIAL DH APPROXIMATED AND THE LEARNT MODEL OF THE REAL MICRO-IGES ROBOT (whose input are the Micro-IGES motor position values), two hidden layers with 30 and 10 neurons each, and one output layer.

TABLE IV RESULTS
FOR THE REAL WORLD CONTROL TESTS.|| P || IS THE MEAN TIP POSITION ERROR NORM (IN MM): IVA) ON THE CIRCULAR PATH; IVB) ON THE SQUARE PATH