Sage Journals: Discover world-class research

Abstract

Minimally invasive surgery (MIS) during cardiovascular interventions reduces trauma and enables the treatment of high-risk patients who were initially denied surgery. However, restricted access, reduced visibility and control of the instrument at the treatment locations limits the performance and capabilities of such interventions during MIS. Therefore, the demand for technology such as steerable sheaths or catheters that assist the clinician during the procedure is increasing. In this study, we present and evaluate a robotically actuated delivery sheath (RADS) capable of autonomously and accurately compensating for beating heart motions by using a model-predictive control (MPC) strategy. We develop kinematic models and present online ultrasound segmentation of the RADS that are integrated with the MPC strategy. As a case study, we use pre-operative ultrasound images from a patient to extract motion profiles of the aortic heart valve (AHV). This allows the MPC strategy to anticipate for AHV motions. Further, mechanical hysteresis in the steering mechanism is compensated for in order to improve tip positioning accuracy. The novel integrated system is capable of controlling the articulating tip of the RADS to assist the clinician during cardiovascular surgery. Experiments demonstrate that the RADS follows the AHV motion with a mean positioning error of 1.68 mm. The presented modelling, imaging and control framework could be adapted and applied to a range of continuum-style robots and catheters for various cardiovascular interventions.

Keywords

Model predictive control robotically actuated delivery sheath ultrasound guided-control beating heart compensation

1. Introduction

According to the World Health Organization, cardiovascular diseases are among the leading causes of death globally (Alwan et al., 2011). Technological developments have the ability to improve and enable treatment of cardiovascular diseases (Himbert et al., 2009; Walther et al., 2009; Wong et al., 2010). Traditionally, open surgery via sternotomy is performed to gain surgical access to the heart, while a heart–lung machine provides life support. Accessibility to the heart is considered to be a major advantage during open surgery, which is at the expense of severe patient trauma. As an alternative, minimally invasive surgery (MIS) could reduce trauma and enable the treatment of high-risk patients who were initially denied surgery (Alfieri et al., 2007; Cribier et al., 2004; Himbert et al., 2009; Seeburger et al., 2010; Ye et al., 2010). However, limited access, vision and control of the instrument at the treatment locations impedes the performance of MIS. During interventions performed without a heart–lung machine, often beating heart motion compensation is desired and could potentially enable future cardiovascular interventions. Tracking of the beating heart requires the attention of the surgeon, whose accuracy deteriorates for complex repetitive motions up to 60 beats per minute (BPM) (Falk, 2002). Therefore, there is a demand for technology that assists the clinician via shared control during the procedure.

In this study, we focus on a model-predictive control (MPC) strategy that anticipates and compensates for the beating heart motion (Figure 1). The MPC strategy is used to compensate the beating heart motion by using a robotically actuated delivery sheath (RADS). Another advantage of MPC is the ability to integrate motion constraints on the instrument in order to prevent damage to sensitive tissue during the procedure. As a case study, we use transcatheter aortic valve implantation (TAVI) to demonstrate the potential of our MPC strategy. Further, applications such as ablation, valve repair surgery and mitral valve chordal replacement can potentially benefit from such an MPC strategy. TAVI can be performed via the transfemoral (TF) or transapical (TA) routes. In this case study, we focus on the TAVI-TA approach, which provides direct surgical access to the aortic and mitral heart valves (Johansson et al., 2011; Ye et al., 2010).

Fig. 1.

Model predictive control (MPC) can be used to steer the robotically actuated delivery sheath (RADS) in order to assist the clinician during cardiovascular surgery. The potential of MPC in cardiovascular surgery can be demonstrated by compensating for the aortic heart valve (AHV) motion in a representative case of transapical transcatheter aortic valve implantation. The RADS ① is inserted through the apex ② into the left ventricle ③ and oriented towards to the aortic annulus ④. The articulating tip ⑤ of the RADS can be steered inside the left ventricle under ultrasound ⑥ image guidance in two degrees of freedom by two pairs of antagonistically configured tension wires. Pre-operative ultrasound data ⑦ can be used as an input to the MPC strategy to anticipate and compensate the AHV motion during surgery.

1.1. Related work

Recent studies have demonstrated the capabilities of MIS in cardiovascular applications such as angioplasty and patent foramen ovale (Gosline et al., 2012; Jayender et al., 2008). These studies describe the use of robotic instruments in order to enable MIS. However, MIS in cardiovascular applications such as coronary artery bypass, ablation and valve surgery could significantly benefit from beating heart compensation. Active robotic stabilization at the treatment location could provide a virtually still scenario that allows the clinician to perform the primary task of the procedure as if the heart was stopped. Compensation of beating heart motions has been extensively covered in research. Studies that focused on the beating heart surface have been described (Bebek and Cavusoglu, 2007; Gangloff et al., 2006; Ginhoux et al., 2005; Nakamura et al., 2001; Ortmaier et al., 2005; Richa et al., 2010; Tuna et al., 2013). These studies cover various predictive strategies by filtering and control to anticipate the beating heart motions. Ortmaier et al. (2005) reported significant correlations between beating heart motions and electrocardiogram (ECG) signals, which were utilized for synchronization by Bebek and Cavusoglu (2007). The majority of these studies used predictive strategies to improve accuracy and to anticipate for occlusions in endoscopic camera images, while this study focuses on ultrasound-guided repair surgery in the beating heart.

Cardiovascular procedures such as ablation and valve-implantation or -repair surgery are often assisted by two-dimensional (2D) or three-dimensional (3D) ultrasound, which are used to evaluate cardiac functionalities (e.g. valve closing) (Walther et al., 2009). The presence of ultrasound imaging during the procedure provides the potential for online soft tissue and instrument tracking. Yuen et al. (2008) evaluated predictive filtering to enable ultrasound-guided robotic tracking of cardiac motions for potential use in ablation. A study by Kesner and Howe (2014) demonstrated capabilities of applying a constant force against a moving target by using ultrasound image feedback and force sensing. Whereas Bowthorpe et al. (2014) focused on predictive control using an ultrasound-guided tele-operated robot to anticipate for cardiac motions. All of these methods describe one dimensional motion compensation, which is along the lateral axis of the instrument. However, existing applications in valve implantation and repair surgery such as mitraclip placement could benefit from motion compensation in two or more degrees of freedom (DOFs) (Tamburino et al., 2010; Wong et al., 2010; Ye et al., 2010). Further, beating heart motion compensation in two or more DOFs could potentially enable future applications in minimally invasive cardiovascular surgery. In this study, we cover an integrated system capable of autonomously and accurately compensating for beating heart motions in two dimensions by using an MPC strategy. Further, our MPC strategy guided by ultrasound imaging considers active constraint handling in order to avoid damage to sensitive tissue. Major aspects in MPC of the integrated system are kinematic modelling of continuum-style robots and feedback on its pose.

Kinematic modelling is used to describe the RADS shape and articulating tip motion. Modelling of continuum-style robots have been investigated by several groups (Bardou et al., 2010; Camarillo et al., 2008; Ding et al., 2010; Dupont et al., 2010; Reiter et al., 2012; Rone and Ben-Tzvi, 2014; Webster and Jones, 2010). Various continuum-style robotic instruments with potential in cardiovascular applications have been described in literature. These devices can be mainly classified into tendon-driven robots and pre-bent concentric tubes. We describe the kinematics of the tendon-driven constant-curvature RADS using both robot-specific and -independent submappings (Webster and Jones, 2010). The kinematics presented in this study considers a rigid link which is attached to the flexible segment of the RADS. However, modelling mismatches and disturbances acting on the system affect the tip positioning performance of continuum-style robots.

Tracking of continuum-style robots plays an important role in order to reduce control tracking errors caused by modelling mismatches and external disturbances. In TAVI, the clinician is often assisted by 2D and 3D transoesophageal or transthoracic echocardiography (Walther et al., 2009). 2D ultrasound has been used for 2D and 3D tracking of flexible instruments by several groups (Hong et al., 2004; Neshat and Patel, 2008; Neubach and Shoham, 2010; Vrooijink et al., 2013, 2014). 3D ultrasound-based tracking of cardiac catheters (some equipped with markers) have also been described in the literature (Koolwal et al., 2010; Nadeau et al., 2015; Novotny et al., 2007; Stoll et al., 2012). In this study, we integrate an online and a robust tracking algorithm that uses 2D ultrasound images of the RADS.

1.2. Contributions

In this study, we demonstrate an integrated system that assists the clinician by compensating for beating heart motions, which could be used in existing and potential future cardiovascular interventions. By active stabilization of the instrument in the beating heart, a virtually still treatment location could be provided. This allows the clinician to perform the procedure as if the heart was stopped. Such a system requires instrument modelling and tracking, which is integrated in a control strategy capable of anticipating for beating heart motions. In this study, we incorporate the forward, inverse and differential kinematic models of the RADS, which considers a rigid link or a tool attached to the flexible tip segment. Subsequently, we integrate a robust 2D ultrasound image segmentation algorithm capable of online evaluation of the RADS tip position. Furthermore, we provide the control objective and corresponding MPC strategy, which considers kinematic modelling and ultrasound feedback. The MPC strategy anticipates for beating heart motions by incorporating models based on pre-operative patient data. In addition to the aforementioned contributions, the MPC strategy considers active constraint handling in order to prevent damage to sensitive tissue. To the best of the authors’ knowledge, we are the first to present an ultrasound-guided MPC strategy capable of compensation for beating heart motions in 2D using a RADS. As a case study, we focus on TAVI-TA by integrating an aortic heart valve (AHV) model in the MPC strategy, which is based on pre-operative 2D ultrasound patient images that describes the AHV motion during the cardiac cycle. The presented strategy for beating heart compensation could be applicable to a wide variety of existing and potential future cardiovascular interventions. Further, we improved the RADS tip positioning accuracy by compensating for instrument hysteresis. Three experimental scenarios using a water container setup demonstrate MPC of the novel integrated system that autonomously and accurately compensates for AHV motions.

2. Methods

This section presents techniques to enable tracking and MPC-based steering of the RADS. Details on the design of the RADS are provided by Vrooijink et al. (2014). We summarize the forward, inverse and differential kinematics which are used in MPC. Furthermore, we describe the 2D ultrasound segmentation method of the RADS used as feedback for MPC. The MPC strategy includes modelling of the heart valve motion based on pre-operative 2D ultrasound data. By combining the kinematic modelling and segmentation of the RADS with MPC, the system is able to anticipate the beating heart motion. Note, in the derivations presented, for notational simplicity, we leave out the discrete time variable denoted by k (except where needed for clarity). Further, the nomenclature describing the variables used in this study is provided in Appendix B.

2.1. Instrument modelling

2.1.1. Forward kinematics

The design and integration of the tendon-driven RADS used in this study is shown in Figure 2. Two pairs of antagonistically configured tension wires are used to actuate the articulating tip of the RADS in two DOFs. Each pair of wires is actuated by a single pulley and controls tip movement in a single DOF (Figure 2). The four-wire design allows for tip movement in two DOFs by using two actuators instead of the three that are required in an instrument with three wires. A kinematic model of the tendon-driven continuum-style robot can be described by a robot-specific and a robot-independent submapping (Webster and Jones, 2010). The robot-specific mapping relates the actuator space to the configuration space. The actuator space is given by the angles of the pulleys ( $ψ_{x}$ and $ψ_{y}$ ), while the configuration space is described by arc parameters such as arc curvature ( $κ$ ), arc plane angle ( $ϕ$ ) and arc length ( $ℓ$ )). The robot-independent mapping transforms the configuration space to the task space (intermediate frame ( $Ψ_{i}$ )).

Fig. 2.

An overview of the robotically actuated delivery sheath (RADS). The articulating tip of the RADS is actuated in two degrees of freedom by two pairs (red and green) of antagonistic-configured tension wires driven by a two pulleys with radii ( $r_{p}$ ) and angles ( $ψ_{x}$ and $ψ_{y}$ ). Three coordinate systems are assigned to describe the tip pose of the RADS: $Ψ_{0}$ is the reference frame fixed to the shaft, $Ψ_{i}$ is the intermediate frame assigned to the arc section and frame ( $Ψ_{t}$ ) is fixed to the articulating tip. Displacement of the tension wires ( $t_{1}, \dots, t_{4}$ ) by $δ_{x}$ and $δ_{y}$ (inset centre) results in instrument bending along the x- and y-axes (frame ( $Ψ_{0}$ )), respectively. The arc of the RADS with parameters, bend angle ( $θ$ ), backbone length ( $ℓ$ ), radius (r) and curvature ( $κ$ ) lies in a plane described by the arc plane (inset right). The orientation of the arc plane about the z-axis of the reference frame ( $Ψ_{0}$ ) is denoted by angle ( $ϕ$ ). Further, the tendon distance to the backbone arc ( $ℓ$ ) is denoted $d_{b}$ . A rigid link (not completely shown) of length ( $l_{t}$ ) is attached to the arc (frame ( $Ψ_{i}$ )) of the RADS.

In order to evaluate the arc parameters of the configuration space, the relation between tendon manipulation at the base and the resulting arc needs to be described. In the derivation presented, we denote for notational simplicity: $c_{*} = \cos (*)$ and $s_{*} = \sin (*)$ . The RADS used in this study is actuated using four tendons ( $t_{i}$ ), where $(i = 1, \dots, 4$ ) with corresponding tendon lengths ( $l_{i}$ ) as shown in Figure 2 (inset centre). The relationship between the arc length of the RADS ( $ℓ$ ) and the arc length of a single tendon ( $l_{i}$ ) is given by

ℓ = l_{i} + θ d_{b} c_{ϕ_{i}}

(1)

where $θ$ is the bend angle, $d_{b}$ denotes the distance between the backbone and tendon and $ϕ_{i}$ describes the angle between the bending direction of the RADS and the location of a single tendon ( $t_{i}$ ). Note that the distance ( $d_{b}$ ) in our instrument is equal for all tendons (Figure 2). Further, the bend angle ( $θ$ ) is related to the curvature by $θ = κ ℓ$ .

Given the tendon configuration as illustrated in Figure 2, we can formulate the individual tendon angles by $c_{ϕ_{1}} = c_{ϕ}$ , $c_{ϕ_{2}} = s_{ϕ}$ , $c_{ϕ_{3}} = - c_{ϕ}$ and $c_{ϕ_{4}} = - s_{ϕ}$ . By combining the individual tendon angles with (1), we can obtain an expression between the arc length ( $ℓ$ ) of the RADS and the individual tendon lengths ( $l_{i}$ ), which yields

ℓ = \frac{l_{1} + l_{2} + l_{3} + l_{4}}{4}

(2)

The antagonistic configuration allows for pairing of the tendons described by $t_{1}$ and $t_{3}$ , and $t_{2}$ and $t_{4}$ . By using the tendon pairs and (2), the arc plane angle ( $ϕ$ ) is obtained according to

ϕ = \arctan (\frac{l_{4} - l_{2}}{l_{3} - l_{1}})

(3)

while the curvature is evaluated as

κ = \frac{(l_{1} - 3 l_{2} + l_{3} + l_{4}) \sqrt{{(l_{4} - l_{2})}^{2} + {(l_{3} - l_{1})}^{2}})}{d_{b} (l_{1} + l_{2} + l_{3} + l_{4}) (l_{4} - l_{2})}

(4)

Note that the arc parameters ( $ℓ$ , $ϕ$ and $κ$ ) described in (2)–(4) are functions of the individual tendon lengths ( $l_{i}$ ). These individual tendon lengths ( $l_{i}$ ) are manipulated by displacements ( $δ_{x}$ ) and ( $δ_{y}$ ), which are introduced by two actuated pulleys (Figure 2). The relation between tendon displacements ( $δ_{x}$ and $δ_{y}$ ) and pulley angles ( $ψ_{x}$ and $ψ_{y}$ ) are given by $δ_{x} = r_{p} ψ_{x}$ and $δ_{y} = r_{p} ψ_{y}$ , where $r_{p}$ describes the pulley radius (equal radii for both pulleys). Hence, we can use (1) to formulate an expression for all tendon lengths ( $l_{i}$ ) as a function of the pulley angles ( $ψ_{x}$ and $ψ_{y}$ ) by $l_{1} = ℓ - r_{p} ψ_{x}$ , $l_{2} = ℓ - r_{p} ψ_{y}$ , $l_{3} = ℓ + r_{p} ψ_{x}$ and $l_{4} = ℓ + r_{p} ψ_{y}$ . Substituting the expressions for the tendon lengths ( $l_{i}$ ) into (3) and (4) yields

ϕ = \arctan (\frac{ψ_{y}}{ψ_{x}})

(5)

and

κ = \frac{r_{p} \sqrt{ψ_{x}^{2} + ψ_{y}^{2}}}{ℓ d_{b}}

(6)

respectively. Hence, we obtain the arc parameters ( $κ$ , $ϕ$ and $ℓ$ ) of the configuration space as a function of the pulley angles ( $ψ_{x}$ and $ψ_{y}$ ) of the actuator space.

The robot-independent mapping is given by the homogeneous transformation matrix ( $H_{i}^{0} \in SE 3$ ), where

H_{i}^{0} = [\begin{matrix} c_{ϕ} c_{κ ℓ} & - s_{ϕ} & c_{ϕ} s_{κ ℓ} & \frac{c_{ϕ} (1 - c_{κ ℓ})}{κ} \\ s_{ϕ} c_{κ ℓ} & c_{ϕ} & s_{ϕ} s_{κ ℓ} & \frac{s_{ϕ} (1 - c_{κ ℓ})}{κ} \\ - s_{κ ℓ} & 0 & c_{κ ℓ} & \frac{s_{κ ℓ}}{κ} \\ 0 & 0 & 0 & 1 \end{matrix}]

(7)

which describes the mapping between the intermediate frame ( $Ψ_{i}$ ) and the reference frame ( $Ψ_{0}$ ).

Remark (Discontinuity for $κ \neq 0$ ). Discontinuities in kinematics can be addressed by substituting differentiable cardinal sine and versine functions according to

sinc (κ ℓ) = {\begin{matrix} \frac{\sin (κ ℓ)}{κ ℓ} & κ ℓ \neq 0 \\ 1 & κ ℓ = 0 \end{matrix}

(8)

and

verc (κ ℓ) = {\begin{matrix} \frac{1 - \cos (κ ℓ)}{κ ℓ} & κ ℓ \neq 0 \\ 0 & κ ℓ = 0 \end{matrix}

(9)

which are defined for $κ = 0$ .

In order to describe the articulating tip frame ( $Ψ_{t}$ ) of the RADS with respect to the intermediate frame ( $Ψ_{i}$ ) as demonstrated in Figure 2, we use a second transformation. This transformation describes the rigid link attached to the arc section (intermediate frame ( $Ψ_{i}$ )) of the RADS, which is given by homogenous transformation matrix ( $H_{t}^{i} \in SE 3$ ), where

H_{t}^{i} = [\begin{matrix} I_{3} & L_{t}^{i} \\ 0_{1 \times 3} & 1 \end{matrix}]

(10)

where $I$ and $0$ represents an identity and a matrix filled with zeros, respectively. Further, $L_{t}^{i} \in R^{3}$ , where $L_{t}^{i} = [0 0 l_{t}]^{T}$ describes the rigid link section by a translation along the z-axes of the intermediate frame ( $Ψ_{i}$ ). We describe the RADS articulating tip pose ( $H_{t}^{0}$ ) in the reference frame ( $Ψ_{0}$ ) by a series of transformations according to $H_{t}^{0} = H_{i}^{0} H_{t}^{i}$ , which completes the forward kinematics.

2.1.2. Inverse kinematics

The inverse kinematics is used to express the pulley angles ( $ψ_{x}$ and $ψ_{y}$ ) as a function of the reference tip position ( $r_{t}^{0} \in R^{3}$ , where $r_{t}^{0} = [r_{x} r_{y} r_{z}]^{T}$ )). In order to evaluate the arc parameters of the configuration-space given a measured or reference tip position ( $r_{t}^{0}$ ), we first determine the arc plane angle ( $ϕ$ ) by

ϕ = \arctan (\frac{r_{y}}{r_{x}})

(11)

Subsequently, we use the forward kinematics presented in (7) and (10) to derive an expression for the reference tip position ( $r_{t}^{0}$ ) according to

[\begin{matrix} r_{t}^{0} \\ 1 \end{matrix}] = H_{t}^{0} [\begin{matrix} o_{t} \\ 1 \end{matrix}] = [\begin{matrix} c_{ϕ} (\frac{(1 - c_{κ ℓ})}{κ} + l_{t} S_{k} ℓ) \\ S_{ϕ} (\frac{(1 - c_{κ ℓ})}{κ} + l_{t} S_{k} ℓ) \\ \frac{S_{k} ℓ}{κ} + l_{t} C_{k} ℓ \\ 1 \end{matrix}]

(12)

where $o_{t} = [0 0 0]^{T} \in R^{3}$ represents the origin of the articulating tip frame ( $Ψ_{t}$ ). Note, that the arc length ( $ℓ$ ) and rigid link length ( $l_{t}$ ) are known. Thus, by substituting (11) into (12), the arc curvature ( $κ$ ) can be solved numerically as described in Algorithm 1.

Algorithm 1 Numerical estimation of the arc curvature ( $κ$ )
Inputs:
$r_{t}^{0} = [r_{x} r_{y} 0]^{T}$ ⊳ Robotically actuated delivery sheath tip position
$ε_{th}$ ⊳ Specified tolerance
Parameters:
$ℓ$ ⊳ Arc length
$l_{t}$ ⊳ Rigid link length
Outputs:
$\hat{κ}$ ⊳ Estimated arc curvature
Method:
1: $ϕ = \arctan (\frac{r_{y}}{r_{x}})$ ⊳ Arc plane angle
2: $\tan (\hat{θ}) = \frac{\| r_{t}^{0} \|}{ℓ + l_{t}}$
3: $\hat{κ} = \frac{2}{ℓ + l_{t}} \tan (\hat{θ})$ ⊳ Initial estimate $\hat{κ}$
4: while $e > ε_{th}$ do
5: $[\begin{matrix} {\hat{r}}_{t}^{0} \\ 1 \end{matrix}] = H_{t}^{0} (\hat{κ}) [\begin{matrix} o_{t} \\ 1 \end{matrix}]$ ⊳ Compute tip position with $\hat{κ}$
6: $e = \| {\hat{r}}_{t}^{0} \| - \| r_{t}^{0} \|$ ⊳ Compute error (x- and y-axes)
7: $\hat{κ} = \hat{κ} + G \cdot e$ ⊳ Re-estimate $\hat{κ}$ , with $G = 4450$ (empirically determined)
8: end while

We first estimate an initial value for the arc curvature ( $\hat{κ}$ ) using trigonometric relations as depicted in Figure 3. Subsequently, we use the initial estimated arc curvature ( $\hat{κ}$ ) to compute the forward kinematics described in (7) and (10). The error (e) between the estimate ( ${\hat{r}}_{t}^{0}$ ) and reference ( $r_{t}^{0}$ ) tip position is evaluated. The computed error (e) is used to correct the estimated arc curvature ( $\hat{κ}$ ), where the feedback gain (G) is included to obtain a desired rate of convergence. The process is repeated until the error (e) is below a specified tolerance ( $ε_{th}$ ). Hence, we obtain an accurate estimation of the arc curvature ( $κ$ ). Given the evaluated arc curvature ( $κ$ ) and arc plane angle ( $ϕ$ ), we can solve (1) for all individual tendon lengths ( $l_{i}$ ) and determine the corresponding pulley angles ( $ψ_{x}$ and $ψ_{y}$ ), which completes the inverse kinematics.

Fig. 3.

Overview of the geometric relations used to estimate an initial value for the arc curvature ( $\hat{κ}$ ) of the robotically actuated delivery sheath (RADS). The arc length ( $ℓ$ ) and rigid link length ( $l_{t}$ ) combined with the RADS tip position ( ${| r}_{t}^{0} |$ ) are used to estimate the bend angle ( $\hat{θ}$ ). The estimated bend angle ( $\hat{θ}$ ) and the approximated length ( $\frac{1}{2} (ℓ + l_{t})$ ) are used to evaluate an initial estimate of the arc curvature ( $\hat{κ}$ ) of the RADS.

2.1.3. Differential kinematics

Differential kinematics is used to relate the change in instrument arc parameters ( $u = [ϕ κ]^{T}$ ) to the change in tip position, which is integrated in the MPC strategy. Similar to the reference tip position in (12), we can describe the tip position according to $p = F (ϕ, κ, ℓ, l_{t})$ , where $p = [p_{x} p_{y} p_{z}]^{T} \in R^{3}$ is the tip position in x-, y- and z-axes (frame ( $Ψ_{0}$ )), which is described by a function ( $F (ϕ, κ, ℓ, l_{t})$ ) of the arc parameters ( $ϕ$ and $κ$ ). Given the instrument design, we consider the arc length ( $ℓ$ ) and rigid link length ( $l_{t}$ ) to be constant in ( $F (ϕ, κ, ℓ, l_{t})$ ). Note, that the arc parameters ( $ϕ$ and $κ$ ) are controlled by pulley angles ( $ψ_{x}$ and $ψ_{y}$ ). An expression that relates the change in instrument arc parameters ( $\overset{\cdot}{u} = [\overset{\cdot}{ϕ} \overset{\cdot}{κ}]^{T}$ ) to instrument tip velocity ( $\overset{\cdot}{p} = [{\overset{\cdot}{p}}_{x} {\overset{\cdot}{p}}_{y} {\overset{\cdot}{p}}_{z}]^{T}$ ) is given by

\overset{\cdot}{p} = J_{R} (u) \overset{\cdot}{u}

(13)

where

\begin{matrix} J_{R} (u) = [\begin{matrix} \frac{\partial F_{1}}{\partial ϕ} & \frac{\partial F_{1}}{\partial κ} \\ \frac{\partial F_{2}}{\partial ϕ} & \frac{\partial F_{2}}{\partial κ} \\ \frac{\partial F_{3}}{\partial ϕ} & \frac{\partial F_{3}}{\partial κ} \end{matrix}] = \\ [\begin{matrix} s_{ϕ} (\frac{c_{κ ℓ} - 1}{κ} - l_{t} s_{κ ℓ}) & c_{ϕ} (\frac{c_{κ ℓ} - 1}{κ^{2}} + \frac{ℓ s_{κ ℓ}}{κ} + ℓ l_{t} c_{κ ℓ}) \\ - c_{ϕ} (\frac{c_{κ ℓ} - 1}{κ} - l_{t} s_{κ ℓ}) & s_{ϕ} (\frac{c_{κ ℓ} - 1}{κ^{2}} + \frac{ℓ s_{κ ℓ}}{κ} + ℓ l_{t} c_{κ ℓ}) \\ 0 & \frac{ℓ c_{κ ℓ}}{κ} - \frac{s_{κ ℓ}}{κ^{2}} - ℓ l_{t} s_{κ ℓ} \end{matrix}] \end{matrix}

(14)

uses the forward kinematics described in (12) to analytically derive the differential kinematics. Note, that similar to (7), discontinuities for $κ = 0$ can be addressed by using derivatives of cardinal sine (sinc $(κ ℓ)$ ) and versine (verc $(κ ℓ)$ ) functions. Given the presented techniques that provide forward, inverse and differential kinematics, often unknown modelling mismatches and external disturbances degrade the control performance. Hence, we include ultrasound feedback in order to improve our control strategy.

2.2. Ultrasound image segmentation

This section elaborates on the segmentation techniques applied to evaluate the RADS centroid location in 2D ultrasound images. In order to view the tip of the RADS in ultrasound images, we insert the instrument in a container and use water as an acoustic transport medium. A radial cross-sectional view of the RADS is obtained by orientating the 2D ultrasound image plane perpendicular to the shaft of the instrument (Figure 4). The RADS tip is positioned in the ultrasound image plane by axial positioning, which is described in details in Section 2.3.4. Note, that by axial positioning, the segmented RADS tip frame in ultrasound images (frame ( $Ψ_{u}$ )) is expressed in the fixed reference frame ( $Ψ_{0}$ ) by position feedback.

Fig. 4.

Ultrasound image segmentation to evaluate the centroid location ( $p_{x}, p_{y}$ ) of the robotically actuated delivery sheath (RADS). The RADS tip (frame ( $Ψ_{t}$ )) is positioned in the ultrasound image plane denoted by frame ( $Ψ_{u}$ ) using axial positioning along the x-axis (frame ( $Ψ_{0}$ )). (a) Radial cross-sectional view of the RADS in 2D ultrasound images. (b) Gaussian filtering using a 2D kernel. (c) Canny edge detection with hysteretic thresholding. (d) Random sample consensus (RANSAC) to localize the centroid ( $p_{x}$ , $p_{y}$ ) of the RADS (centre of the blue circle). The green and red points are considered inliers and outliers, respectively.

A representative ultrasound image of the instrument tip is shown in Figure 4(a). We observe a semi-circular shape that describes the reflecting surface of the RADS. By evaluation of the pre-operative ultrasound data (Figure 6), we did not observe identical semi-circular shapes that represented anatomical structures. Therefore, we use segmentation techniques based on the circular shape parameters to localize the tip of the RADS in ultrasound images. Further, the acoustic impedance difference between RADS materials and cardiac tissue is considered to be significantly large, which enhances the visibility of the instrument (Aldrich, 2007; Vrooijink et al., 2014). This would enable instrument detection when interaction with cardiac tissue is considered.

Preprocessing is performed to enhance the visibility of the RADS in ultrasound images. In order to reduce speckle and to smoothen edges in the ultrasound image, a 2D Gaussian kernel is applied, as depicted in Figure 4(b). The contrast between instrument and the environment in ultrasound images is sufficient for edge detection. Hence, we use a Canny edge detector with hysteresis thresholding to obtain an edge map of the ultrasound image as shown in Figure 4(c) (Forsyth and Ponce, 2002). Hysteretic thresholding reduces the detection of irrelevant edges, that do not describe the semi-circular surface of the RADS. However, irrelevant edges, often introduced by surface deformations caused by artifacts and bending of the instrument may still exist. Hence, the centroid location is evaluated using a random sample consensus (RANSAC) strategy (Algorithm 2). RANSAC is robust to irrelevant edges that do not fit the parametric description of the semi-circular model representing the RADS (Forsyth and Ponce, 2002).

Algorithm 2 Random sample consensus robotically actuated delivery sheath localization
Inputs:
$A_{c} \leftarrow {x_{c, v \| v = 1, \dots, w}}$ ⊳ Set of detected edge points ( $x_{c}$ )
$f : H_{3} \to m_{c}$ ⊳ Computes the algebraic circle model parameters ( $m_{c}$ ) from a set ( $H_{3}$ ) of three randomly selected edge points
$C (m_{c}, x_{c})$ ⊳ Cost function for a single edge point (1 if $x_{c}$ is an inlier to the algebraic circle parameters ( $m_{c}$ ), 0 otherwise)
n ⊳ Number of iterations
Outputs:
$m_{c}^{*}$ ⊳ Best model parameters
$S_{c}^{*}$ ⊳ Best consensus set (inliers)
$J_{c}^{*}$ ⊳ Best cost
Method:
1: for $i \leftarrow 1, n$ do
2: $H_{3, i} \leftarrow random_3 pnts (A_{c})$ ⊳ (I) Hypothesis
3: $m_{c}, i \leftarrow f (H_{3, i})$
4: if $suffice (m_{c, i})$ then ⊳ (II) Preliminary test
5: $S_{c}, i \leftarrow {\forall x_{c} \in A_{c} \| C (m_{c, i}, x_{c}) = 1}$
6: $J_{c}, i \leftarrow \sum_{x_{c} \in A_{c}} C (m_{c, i}, x_{c})$
7: if $J_{c}^{*} < J_{c}, i$ then ⊳ (III) Evaluation
8: $J_{c}^{*} \leftarrow J_{c, i}$
9: $m_{c}^{*} \leftarrow m_{c, i}$
10: $S_{c}^{*} \leftarrow S_{c, i}$
11: end if
12: end if
13: $i \leftarrow i + 1$
14: end for

The evaluated set ( $A_{c}$ ) of edge points ( $x_{c} \in R^{2}$ ) from the Canny edge detector is used as an input to the RANSAC algorithm. The RANSAC algorithm is an iterative process that consists of three steps (i.e. (I) hypothesis, (II) preliminary test and (III) evaluation). In (I), a hypothetical solution is found by fitting algebraic circle model parameters (m) for $f : H_{3} \to m_{c}$ to a set ( $H_{3}$ ) of candidate inliers. The set ( $H_{3}$ ) consists of three arbitrary selected edge points from the evaluated set ( $A_{c}$ ). In (II), a preliminary test $(suffice (m_{c, i}))$ is performed to evaluate the model parameters. The model parameters such as radius and centroid location should be consistent with those of the RADS in order to qualify as a potential solution. Further, all edge points of the edge map are evaluated using a cost function $(C (m_{c}, x_{c}))$ to determine whether they are sufficiently close to the periphery of the circular shape. The fitted model is acceptable if a sufficiently large part of the semi-circular surface has been evaluated as the consensus set. In (III), an evaluation is performed in order to improve the previous solution. The model parameters and consensus set are both refined if the computed cost ( $J_{c, i}$ ) of the current iteration exceeds the previous solution. After n iterations, the centroid location ( $p_{x}, p_{y}$ ) of RADS is evaluated from the model parameters ( $m_{c}$ ) and displayed as the centre of the circle (Figure 4(d)).

The performance of the RANSAC algorithm is evaluated and validated experimentally using a sequence of 600 ultrasound images. Experiments show that the localization error of the RADS decreases if the number (n) of iterations of the RANSAC algorithm are increased as depicted in Figure 5(b). In order to minimize computational costs and considering the ultrasound imaging resolution of approximately 0.12 mm per pixel, the number of RANSAC iterations should be limited to approximately 800, which is depicted green in Figure 4. Further, we compare the results of RADS segmentation for varying RANSAC iterations using manual segmentation as a ground truth (Figure 5(a)). The results show, that by increasing the RANSAC iterations, the segmentation error remains approximately constant (0.4 mm). Since there is no reason to believe that an increase in RANSAC iterations would deteriorate segmentation accuracy, we estimate that the insignificant increased error in Figure 5(a), could be attributed to imperfections in manual segmentations.

Fig. 5.

Experimental evaluation and validation of image segmentation using a sequence of 600 ultrasound images. (a) Compares the performance between varying random sample consensus (RANSAC) iterations using a ground truth obtained by manual segmentation. (b) Describes the relation between the number of iterations and the localization error using a ground truth obtained after $10^{6}$ RANSAC iterations.

The presented ultrasound segmentation method is used in our closed-loop control strategy to improve the performance. Note, that the presented 2D segmentation method can be expanded to 3D ultrasound imaging by parallel evaluation of slices along the instrument shaft.

2.3. Model predictive control

In this section we present the MPC architecture used to control the RADS in closed loop. The objective of the MPC architecture is to anticipate the AHV motion based on a model. This allows the clinician to focus on the primary task of the procedure, while compensation for AHV motions is provided. In order to develop an AHV model, we use the pre-operative 2D ultrasound images that describe the AHV motion of a human during the cardiac cycle. Further, we integrate forward, inverse and differential kinematics described in Section 2.1 to determine instrument motion. The 2D ultrasound segmentation presented in Section 2.2 is used to provide instrument feedback for control.

2.3.1. Model description

The continuous time representation of the model used in the MPC strategy incorporates the differential kinematics described in (13). A cascade configuration is used according to

{\overset{\cdot}{p}}_{(t)} = J_{R} (u_{(t)}) {\overset{\cdot}{u}}_{(t)}

(15)

{\overset{\cdot\cdot}{u}}_{(t)} = v_{(t)} + e_{(t)}

(16)

y_{(t)} = p_{(t)} + e_{(t)}

(17)

where $v_{(t)}$ is considered to be the MPC control input signal. However, the change in arc parameters ( $\overset{\cdot}{u_{(t)}}$ ) is used to steer the RADS according to the MPC strategy. The system is subject to an inequality constraint according to

p_{\min} \leq p_{(t)} \leq p_{\max} for all t

(18)

where the instrument tip position is restricted. The instrument tip position is limited in order to prevent damage to the instrument arc section. Further, by adding restricted regions for the instrument tip position, damage to sensitive tissue can be avoided. In order to preserve dominant features required for ultrasound image segmentation, we add an additional inequality constraint. The constraint is used to limit the instrument tip velocity in x- and y-directions of frame ( $Ψ_{0}$ ) by

{\overset{\cdot}{p}}_{\min} \leq {\overset{\cdot}{p}}_{(t)} \leq {\overset{\cdot}{p}}_{\max} for all t

(19)

The discrete time model integrated within the MPC strategy can be obtained by using forward Euler discretization of the continuous time model according to

p_{(k + 1)} = p_{(k)} + J_{R} (u_{(k)}) a_{(k)}

(20)

a_{(k + 1)} = a_{(k)} + T_{s}^{2} v_{(k)} + T_{s}^{2} e_{(k)}

(21)

y_{(k)} = p_{(k)} + e_{(k)}

(22)

where we substitute $a_{(k)} = u_{(k + 1)} - u_{(k)}$ and $T_{s}$ denotes the sampling time. Similarly, the equivalent discrete time constraints are given by

p_{\min} \leq p_{(k)} \leq p_{\max} for all k

(23)

and

{\overset{\cdot}{p}}_{\min} \leq \frac{p_{(k + 1)} - p_{(k)}}{T_{s}} \leq {\overset{\cdot}{p}}_{\max} for all k .

(24)

Note, that the constraints stacked into a vector can be evaluated component-wise. The presented discrete time model is integrated within MPC strategy to optimize the control objective.

2.3.2. Control objective

The main MPC objective is to anticipate the AHV motion, while considering a desired reference input. In order to capture these requirements, we can formulate a reference tracking objective according to

r = r_{AHV} + r_{d}

(25)

where $r_{d}$ describes the RADS desired tip position without a priori knowledge, while $r_{AHV}$ represents the position of the AHV with a priori knowledge obtained from pre-operative modelling. In this study, we assume that a priori knowledge on AHV position is sufficiently well-described by modelling, which is presented in Section 2.3.3. In order to achieve our control objective in MPC, we formulate a generalized predictive control (GPC) cost function (J) (van den Boom and Stoorvogel, 2012). This cost function evaluates, over a horizon, the reference tracking error and the control action by

\begin{matrix} J_{(v, k)} = \sum_{j = N_{m}}^{N} {({\hat{y}}_{(k + j | k)} - r_{(k + j)})}^{T} ({\hat{y}}_{(k + j | k)} - r_{(k + j)}) \\ + \sum_{j = 1}^{N} v_{(k + j - 1 | k)}^{T} λ^{T} λ v_{(k + j - 1 | k)} \end{matrix}

(26)

where $N = 15$ describes the prediction horizon, $N_{m}$ is the minimum cost horizon and equals one, ( ${\hat{y}}_{(k + j | k)}$ ) is the predicted instrument tip position ( $y_{(k + j)}$ ) based on knowledge up to time (k) and $λ$ is the control input weighting matrix given by

λ = [\begin{matrix} 1 \times 10^{- 2} & 0 \\ 0 & 1 \times 10^{- 7} \end{matrix}]

(27)

which evaluates the control action ( $v$ ). The control input weighting matrix ( $λ$ ), the prediction horizon (N) and minimum cost-horizon ( $N_{m}$ ) are empirically determined. By rewriting the cost function (26) as

J_{(v, k)} = \sum_{j = 0}^{N - 1} {\hat{z}}_{(k + j | k)}^{T} Γ_{(j)} {\hat{z}}_{(k + j | k)}

(28)

we obtain the standard form used in MPC, where

z_{(k)} = [\begin{matrix} {\hat{y}}_{(k + 1 | k)} - r_{(k + 1)} \\ λ v_{(k)} \end{matrix}]

(29)

describes the reference tracking error and control action, while

Γ_{(j)} = {\begin{matrix} [\begin{matrix} 0_{3 \times 3} & 0_{3 \times 2} \\ 0_{2 \times 3} & I_{2} \end{matrix}] & for 0 \leq j < N_{m} - 1 \\ [\begin{matrix} I_{3} & 0_{3 \times 2} \\ 0_{2 \times 3} & I_{2} \end{matrix}] & for N_{m} - 1 \leq j < N - 1 \end{matrix}

(30)

is a diagonal selection matrix ( $Γ_{(j)} \in R^{5 \times 5}$ ) used to describe the horizon of the cost function. By minimizing the cost function ( $J_{(v, k)}$ ), we optimize between RADS tracking accuracy and control effort, where the reference ( $r_{AHV}$ ) is given by AHV modelling.

2.3.3. AHV modelling

Modelling is used to anticipate the AHV motion in MPC. As a case study, we observe the motion of the AHV annulus during the human cardiac cycle. We use pre-operative 2D transoesophageal echocardiogram (TEE) to obtain the motion profile of a human AHV annulus (Figure 6). The motion profile obtained from 2D ultrasound images is used to demonstrate compensation of the AHV motion. In this study, we assume a constant shape of the aortic annulus. Hence, we obtain the centre position of the aortic annulus by manual segmentation (Figure 6).

Fig. 6.

A pre-operative two-dimensional transoesophageal echocardiogram (TEE) of the aortic heart valve (AHV) annular plane. The centre location of the aortic heart valve annulus is manually segmented for multiple cardiac cycles. The human evaluated annulus centre positions are used for fitting a motion model. The blue line demonstrates the manually segmented annulus position during the cardiac cycle. The corresponding red line describes the position of the fitted model. Further, the velocities and accelerations of the AHV models are provided.

The manual segmented aortic annulus position during the cardiac cycle is used to derive a model capable of describing the periodic motion. A two-term Fourier series according to

f_{x} = a_{0_{x}} + a_{1_{x}} c_{\frac{k}{S_{r}} ω} + b_{1_{x}} s_{\frac{k}{S_{r}} ω} + a_{2_{x}} c_{2 \frac{k}{S_{r}} ω} + b_{2_{x}} s_{2 \frac{k}{S_{r}} ω}

(31)

and

f_{y} = a_{0_{y}} + a_{1_{y}} c_{\frac{k}{S_{r}} ω} + b_{1_{y}} s_{\frac{k}{S_{r}} ω} + a_{2_{y}} c_{2 \frac{k}{S_{r}} ω} + b_{2_{y}} s_{2 \frac{k}{S_{r}} ω}

(32)

is used to describe the periodic aortic annulus motion, where k describes discrete time, $S_{r}$ denotes the number of samples per second and $ω$ represents the frequency of the periodic function in radians per second. The frequency ( $ω$ ) is given by $ω = \frac{ω_{x} + ω_{y}}{2}$ , and relates to the heart rate ( $H_{r}$ ) in beats per minute (BPM) according to $H_{r} = \frac{60 ω}{2 π}$ . The corresponding coefficients and frequencies ( $a_{0_{*}}$ , $a_{1_{*}}$ , $b_{1_{*}}$ , $a_{2_{*}}$ , $b_{2_{*}}$ and $ω_{*}$ ) are provided in Table 1. By considering the frequencies provided in Table 1, the observed heart rate ( $H_{r}$ ) is 42 BPM. The goodness of data fitting is given by the coefficients of determination $R^{2} = 0.86$ and $R^{2} = 0.97$ along the x- and y-axes, respectively. Note, that an increase in the number of terms of the Fourier series described in (31) and (32) does not result in a significant improvement in the coefficients of determination. From the manual segmentation we observe displacements of $9.60 \pm 3.23$ mm in x-axis, while the y-axis shows displacements of $12.25 \pm 1.27$ mm. Further, from the fitted model we observe maximum velocities and accelerations of approximately 40 mm/s and 250 mm/ $s^{2}$ , respectively. The periodic aortic annulus motions described by $f_{x}$ and $f_{y}$ provides the tracking reference ( $r_{AHV} = [f_{x} f_{y} 0]^{T}$ ) for MPC. In order to anticipate the motion of the aortic annulus, a priori knowledge of the tracking reference ( $r_{AHV}$ ) is considered during experimental evaluation using MPC.

Table 1.

Fourier coefficients used in (31) and (32) to describe the periodic aortic heart valve motion depicted in Figure 6. Subscript * denotes the corresponding x- or y-axes.

	$a_{0_{*}}$	$a_{1_{*}}$	$b_{1_{*}}$	$a_{2_{*}}$	$b_{2_{*}}$	$ω_{*}$
x-axis	0.0547	−1.7850	−2.4530	−0.5737	−0.7064	4.4050
y-axis	0.3157	−3.7290	−2.7090	−2.4190	0.4965	4.4240

2.3.4. Controller design

The MPC strategy used to control the articulating tip of the RADS is presented in Figure 7. The control variable ( $v$ ) is used to determine the change in arc parameters ( $a$ ), which are implemented in order to position the tip of the RADS. Subsequently, the tip position ( $y = [p_{x} p_{y} 0]^{T}$ ) of the RADS is measured using a 2D ultrasound transducer as described in Section 2.2. However, ultrasound images are often prone to noise and the tip of the RADS may not always be detected during tracking. Hence, we add an extended Kalman filter to provide state estimation ( $\hat{y} = [{\hat{p}}_{x} {\hat{p}}_{y} {\hat{p}}_{z}]^{T}$ ) based on kinematics described in Section 2.1 and 2D ultrasound measurements (Bar-Shalom et al., 2001). The position error ( $e_{rr} = y - \hat{y}$ ) is evaluated in order to adapt the Kalman filter. Note, that the articulating tip of the RADS must intersect with the ultrasound image plane in order to provide position feedback. Hence, we limit the articulating tip motion to the ultrasound image plane by autonomous axial positioning of the RADS along the z-axis of frame ( $Ψ_{0}$ ) using a linear stage. The state estimate ( ${\hat{p}}_{z}$ ) provided by the Kalman filter is used as an input to axial positioning.

Fig. 7.

Model-predictive control (MPC) strategy is used to steer the robotically actuated delivery sheath (RADS). The referenced tip position by the clinician is denoted by $r_{d}$ , while $r_{AHV}$ describes the position of the aortic heart valve with a priori knowledge obtained from pre-operative ultrasound images. The tip position obtained by ultrasound image segmentation is denoted by $y$ , with corresponding filtered position described by $\hat{y}$ . The arc parameters denoted by $a$ is provided as an input to the RADS and the Kalman filter. The resulting positioning error is given by $e_{rr}$ , which is used to adapt the Kalman filter.

In order to describe the MPC strategy, we formulate a constrained standard predictive control problem (CSPCP), which uses a multiple-input and a multiple-output (MIMO) state-space representation of the system (Kinnaert, 1989; van den Boom and Stoorvogel, 2012). We can write the CSPCP using a MIMO model of the RADS in state-space realization according to the following:

x_{(k + 1)} = A x_{(k)} + B_{1} e_{(k)} + B_{2} w_{(k)} + B_{3} v_{(k)}

(33)

y_{(k)} = C_{1} x_{(k)} + D_{11} e_{(k)} + D_{12} w_{(k)} + D_{13} v_{(k)}

(34)

z_{(k)} = C_{2} x_{(k)} + D_{21} e_{(k)} + D_{22} w_{(k)} + D_{23} v_{(k)}

(35)

χ_{(k)} = C_{4} x_{(k)} + D_{41} e_{(k)} + D_{42} w_{(k)} + D_{43} v_{(k)} \leq X_{(k)}

(36)

where we use the model described in (20) to (22). The state is given by $x_{(k)} = {[\begin{matrix} p_{(k)} & a_{(k)} \end{matrix}]}^{T}$ , $e_{(k)}$ represents zero mean white noise, while $w_{(k)} = {[\begin{matrix} d_{(k)} & r_{(k + 1)} \end{matrix}]}^{T}$ combines all known external signals such as disturbance ( $d_{(k)}$ ) and reference ( $r_{(k + 1)}$ ). Further, the generalized input is described by $v_{(k)}$ . The state-space matrices in (33) and (34) are given by

\begin{matrix} A = [\begin{matrix} I_{3} & J_{R} (u) \\ 0_{2 \times 3} & I_{2} \end{matrix}], B_{1} = [\begin{matrix} 0_{3 \times 3} & 0_{3 \times 2} \\ 0_{2 \times 3} & T_{s}^{2} I_{2} \end{matrix}], \\ B_{2} = [\begin{matrix} T_{s} I_{3} & 0_{3 \times 3} \\ 0_{2 \times 3} & 0_{2 \times 3} \end{matrix}], B_{3} = [\begin{matrix} 0_{3 \times 2} \\ T_{s}^{2} I_{2} \end{matrix}], C_{1} = [\begin{matrix} C_{0} 0_{3 \times 2} \end{matrix}], \\ D_{11} = [\begin{matrix} I_{2} 0_{3 \times 2} \end{matrix}], D_{12} = 0_{3 \times 6}, D_{13} = 0_{3 \times 2} \end{matrix}

(37)

where we use the state ( $p$ ) and the inverse kinematics described in Section 2.1.2 to provide the arc parameters ( $u$ ) for the differential kinematics incorporated in matrix ( $A$ ). Further, the output of the system is limited to the 2D ultrasound image plane, which is given in x- and y directions of frame ( $Ψ_{0}$ ) by

C_{0} = [\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{matrix}]

(38)

In this study, we do not model external disturbances, hence we denote $d_{(k)} = [0 0 0]^{T}$ , for all k. The matrices $C_{2}$ , $D_{21}$ , $D_{22}$ and $D_{23}$ associated with the cost signal ( $z_{(k)}$ ) introduced in (35) can be obtained by substituting ${\hat{y}}_{(k + 1 | k)} = C_{1} x_{(k + 1)} + D_{11} {\hat{e}}_{(k + 1 | k)}$ and $v_{(k)}$ into (29), which can be rewritten as

\begin{array}{l} {\hat{z}}_{(k)} = \underset{C_{2}}{\underset{︸}{[\begin{matrix} C_{1} A \\ 0_{2 \times 6} \end{matrix}]}} x_{(k)} + \underset{D_{21}}{\underset{︸}{[\begin{matrix} C_{1} B_{1} \\ 0_{2 \times 3} \end{matrix}]}} e_{(k)} \\ + \underset{D_{22}}{\underset{︸}{([\begin{matrix} C_{1} B_{2} \\ 0_{2 \times 6} \end{matrix}] + [\begin{matrix} 0_{3 \times 3} & - I_{3 \times 3} \\ 0_{2 \times 3} & 0_{2 \times 3} \end{matrix}])}} w_{(k)} + \underset{D_{23}}{\underset{︸}{[\begin{matrix} C_{1} B_{3} \\ λ \end{matrix}]}} v_{(k)} \end{array}

(39)

where we use the zero mean white noise estimate ( ${\hat{e}}_{(k + 1 | k)} = 0_{3 \times 1}$ ).

The constraints presented in (23) and (24) can be combined and described by two one-sided constraint signals ( $χ_{1, (k)}$ and $χ_{2, (k)}$ ) stacked together according to

\begin{array}{l} \underset{χ_{(k)}}{\underset{︸}{[\begin{matrix} χ_{1, (k)} \\ χ_{2, (k)} \end{matrix}]}} = [\begin{matrix} x_{(k + 1)} \\ - x_{(k + 1)} \end{matrix}] = \underset{C_{4}}{\underset{︸}{[\begin{matrix} A \\ - A \end{matrix}]}} x_{(k)} + \underset{D_{41}}{\underset{︸}{[\begin{matrix} B_{1} \\ - B_{1} \end{matrix}]}} e_{(k)} + \\ \underset{D_{42}}{\underset{︸}{[\begin{matrix} B_{2} \\ - B_{2} \end{matrix}]}} w_{(k)} + \underset{D_{43}}{\underset{︸}{[\begin{matrix} B_{3} \\ - B_{3} \end{matrix}]}} v_{(k)} \leq \underset{X_{(k)}}{\underset{︸}{[\begin{matrix} X_{1, m a x} \\ - X_{2, m i n} \end{matrix}]}} \end{array}

(40)

where $X_{1, \max} = {[\begin{matrix} p_{\max} & T_{s} J_{R}^{+} (u) {\overset{\cdot}{p}}_{\max} \end{matrix}]}^{T}$ and $X_{2, \min} = {[\begin{matrix} p_{\min} & T_{s} J_{R}^{+} (u) {\overset{\cdot}{p}}_{\min} \end{matrix}]}^{T}$ . Note, that by using the differential kinematic relation in (13), the constraint on the instrument tip velocity described in (24) can be rewritten as a constraint on the change in arc parameters ( $a_{(k)}$ ) according to

T_{s} J_{R}^{+} (u) {\overset{\cdot}{p}}_{\min} \leq a_{(k)} \leq T_{s} J_{R}^{+} (u) {\overset{\cdot}{p}}_{\max} for all k

(41)

where $J_{R}^{+} (u)$ is the Moore–Penrose pseudo-inverse of $J_{R} (u)$ . By adding the stacked inequality constraint, we complete the CSPCP presented in (33) to (36).

The CSPCP is used to formulate a model based on the concept of prediction, which is described in Appendix C. The prediction model is used to solve for the optimal control vector ( ${\tilde{v}}_{(k)}$ ). By using predictions of the cost signal provided in (39) and the diagonal selection matrix described in (30), given the prediction interval $0 \leq j \leq N - 1$ , we can formulate a signal vector ( ${\tilde{z}}_{k}$ ) similar to (52) and provide a block diagonal selection matrix ( $\bar{Γ}$ ) by

{\tilde{z}}_{(k)} = [\begin{matrix} {\hat{z}}_{(k | k)} \\ {\hat{z}}_{(k + 1 | k)} \\ ⋮ \\ {\hat{z}}_{(k + N - 1 | k)} \end{matrix}] and \bar{Γ} = [\begin{matrix} Γ_{(1)} & 0_{5 \times 5} & \dots & 0_{5 \times 5} \\ 0_{5 \times 5} & Γ_{(2)} & ⋱ & ⋮ \\ ⋮ & ⋱ & ⋱ & 0_{5 \times 5} \\ 0_{5 \times 5} & \dots & 0_{5 \times 5} & Γ_{(N - 1)} \end{matrix}]

(42)

respectively, which can be used to rewrite the cost-function presented in (28) according to

J_{(v, k)} = {\tilde{z}}_{(k)}^{T} \bar{Γ} {\tilde{z}}_{(k)}

(43)

The predicted cost signal ( ${\tilde{z}}_{k}$ ) in (42) can be formulated as

{\tilde{z}}_{(k)} = {\tilde{C}}_{2} x_{(k)} + {\tilde{D}}_{21} e_{(k)} + {\tilde{D}}_{22} {\tilde{w}}_{(k)} + {\tilde{D}}_{23} {\tilde{v}}_{(k)}

(44)

where matrices ${\tilde{C}}_{2}$ , ${\tilde{D}}_{21}$ , ${\tilde{D}}_{22}$ and ${\tilde{D}}_{23}$ can be obtained according (54) and (55). In addition to the constraint described in CSPCP according to (40), we add additional constraints to the prediction model to shape the control signal. In order to obtain a smooth and robust control action, we add a control horizon constraint to the prediction model. The equality constraint ( $ϒ_{(k)} = 0_{2 \times 1}$ ) on the control horizon can be described by

v_{(k + j | k)} = ϒ_{(k)}, for N_{c} \leq j < N

(45)

where the control action is assumed to be zero after control horizon ( $N_{c} = 10$ ), which is empirically determined. The corresponding prediction signal ( ${\tilde{υ}}_{(k)}$ ) according to (51) and (52) is given by

{\tilde{υ}}_{(k)} = {\tilde{C}}_{3} x_{(k)} + {\tilde{D}}_{31} e_{(k)} + {\tilde{D}}_{32} {\tilde{w}}_{(k)} + {\tilde{D}}_{33} {\tilde{v}}_{(k)} = {\tilde{ϒ}}_{(k)}

(46)

where the equality constraint prediction vector is given by ${\tilde{ϒ}}_{(k)} = 0_{10 \times 1}$ and the corresponding matrices are described by ${\tilde{C}}_{3} = 0_{10 \times 6}$ , ${\tilde{D}}_{31} = 0_{10 \times 3}$ , ${\tilde{D}}_{32} = 0_{10 \times 90}$ and

{\tilde{D}}_{33} = [\begin{matrix} 0_{2 \times 2} & \dots & 0_{2 \times 2} & I_{2} & 0_{2 \times 2} & \dots & 0_{2 \times 2} \\ 0_{2 \times 2} & \dots & 0_{2 \times 2} & 0_{2 \times 2} & I_{2} & ⋱ & ⋮ \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋱ & 0_{2 \times 2} \\ 0_{2 \times 2} & \dots & 0_{2 \times 2} & 0_{2 \times 2} & \dots & 0_{2 \times 2} & I_{2} \end{matrix}]

(47)

Further, a prediction vector ( ${\tilde{χ}}_{(k)}$ ) of the constraint signal ( $χ_{(k)}$ ) described in (40) can be constructed similar to (52), with corresponding prediction matrices computed according to (54) and (55). The prediction of the inequality constraint is given by

{\tilde{χ}}_{(k)} = {\tilde{C}}_{4} x_{(k)} + {\tilde{D}}_{41} e_{(k)} + {\tilde{D}}_{42} {\tilde{w}}_{(k)} + {\tilde{D}}_{43} {\tilde{v}}_{(k)} \leq {\tilde{X}}_{(k)}

(48)

which also completes the prediction model.

By minimizing the following cost function

min_{\tilde{v}} {\tilde{z}}_{(k)}^{T} \bar{Γ} {\tilde{z}}_{(k)}

(49)

subject to constraints described in (46) and (48), we can solve for the optimal control vector ( ${\tilde{v}}_{(k)}$ ) that optimizes the CSPCP. The optimization described in (49) can be evaluated as a quadratic programming problem subject to constraints. For details on the derivations, we refer the reader to work of van den Boom and Stoorvogel (2012). We use qpOases software in C++ to solve the quadratic programming problem online (Ferreau et al., 2014). The solution is used to obtain the optimal control vector ( ${\tilde{v}}_{(k)}$ ). We implement according to the receding horizon principle, where we apply ( $a_{(k)}$ ) to steer the RADS according to the optimization.

2.3.5. Hysteresis compensation

Before the optimal solution is implemented, we compensate for hysteresis in the system. Hysteresis caused by backlash often occurs in cable-driven instruments such as endoscopes and catheters (Reilink et al., 2013). Hysteresis can be described by a positive or negative contact mode and a free mode. In contact mode, kinematics describes the relation between pulley angles ( $ψ_{x}$ and $ψ_{y}$ ) and RADS tip position ( $p$ ), which is not known in free mode. Further, the positive and negative contact modes corresponds to the direction in which the pulleys are engaged. Hysteresis often introduces inaccuracies in the response of the instrument. Hence, we use a classical approach to reduce hysteresis in the system according to

Δ c = {[\begin{matrix} Δ ψ_{x} & Δ ψ_{y} \end{matrix}]}^{T} + {[\begin{matrix} f_{(Δ ψ_{x})} & f_{(Δ ψ_{y})} \end{matrix}]}^{T}

(50)

where $Δ c$ is the compensated angular velocity of the pulleys ( $Δ ψ_{x}$ and $Δ ψ_{y}$ ) and

f_{(Δ ψ_{* (k)})} = {\begin{matrix} sign (Δ ψ_{* (k)}) (α_{*}^{+} - α_{*}^{-}) & If sign (Δ ψ_{* (k)}) \neq sign (Δ ψ_{* (k - 1)}) \\ 0 & Otherwise \end{matrix}

is the function that represents the compensation term. In order to limit chatter, we only apply the compensation term if the velocity changes direction and exceeds a threshold. Compensation is provided by hysteretic parameters ( $α_{*}^{+}$ and $α_{*}^{-}$ ), which are the positive and negative contact positions, respectively. The positive and negative contact positions are experimentally evaluated by measuring the hysteresis between pulley angles and corresponding tip locations.

3. Experiments

In this section, we present experiments performed to evaluate the MPC strategy. First, we describe the components and layout of the experimental setup used to steer the RADS. Subsequently, a pre-operative AHV model is presented and integrated with the MPC strategy. This is followed by the experimental plan and results that demonstrate the capabilities of the MPC strategy.

3.1. Experimental setup and materials

The experimental setup used to evaluate the performance of the MPC is shown in Figure 8. The RADS used in our experiments is based on a cable-ring structure surrounded by a hinged tube (DEAM Corporation, Amsterdam, The Netherlands) (Breedveld, 2010). Two pairs of antagonistically configured tension wires facilitate articulating tip movement of the RADS in two DOFs (Vrooijink et al., 2014). The tension wires of the RADS are driven by an ECMax22 motor with a GP32/22 gearhead (Maxon Motor, Sachseln, Switzerland). Further, an LX30 translational stage (Misumi Group Inc., Tokyo, Japan) is used to translate the RADS along the longitudinal axis in order to compensate for the tip motion in the out-of-image-plane direction. The maximum velocity and acceleration along the longitudinal axis are 230 mm/s and 2800 mm/s², respectively. All motors are actuated by an Elmo Whistle 2.5/60 motor controller (Elmo Motion Control Ltd, Petach-Tikva, Israel). The RADS is inserted in a container filled with water in order to enable ultrasound-based tip tracking. A Siemens 18L6 transducer operating with a frequency of 16 MHz, a power level of −4 dB and a scanning depth of 4 cm on a Siemens Acuson S2000 ultrasound system (Siemens AG, Erlangen, Germany) is positioned at the tip of the RADS (in the container) to obtain ultrasound images as feedback. Note, that without any modifications to ultrasound image segmentation, the transthoracic echocardiography approach used in this study could be replaced by TEE. The ultrasound images are transferred via S-video output (in-plane resolution of approximately 0.12 mm per pixel) to a computer (2.80 GHz Intel i7, 8 GB internal memory and 64-bit Windows 7) with a frame rate of 25 frames per second. Further, this computer is used to implement and execute in C++ the MPC strategy which provides control signals to the motors and electronics. The MPC strategy uses a sampling time ( $T_{s}$ ) of 0.04 s. In order to anticipate the beating heart motion, we integrate an AHV model within the MPC strategy.

Fig. 8.

The experimental setup used to control the robotically actuated delivery sheath (RADS) using a model-predictive control strategy. ① RADS. ② Container filled with water in which the RADS is inserted. ③ Ultrasound transducer. ④ Ultrasound image with a radial cross-sectional view of the RADS. ⑤ Motors and corresponding electronics used to control the articulating tip of the RADS. ⑥ Translation along the longitudinal axis of the RADS in order to position the tip in the two-dimensional ultrasound image plane. The top inset shows the flexible segment (articulating tip) of the RADS, which uses a hinged tube construction. The bottom inset depicts a longitudinal cross-section with dimensions of the RADS. An antagonistic configuration of a pair of tension wires (red) is actuated by a pulley-driven system. Each pair of tension wires (total of two pairs) is guided through the flexible shaft and through two incompressible brass tubes (yellow) to actuate a single degree of freedom of the articulating tip.

3.2. Experimental scenarios

A series of experiments have been conducted in order to evaluate the tracking accuracy of the MPC strategy in an integrated system using a reference signal described in (25). In previous research, a model-based approach was used to steer the RADS along a circular path using 2D ultrasound images as feedback (Vrooijink et al., 2014). The results of these experiments showed, while moving at 2 mm/s, mean positioning errors of approximately 2 mm along the x- and y-axes. In this study, we aim to improve these results and provide novel functionalities such as compensation for beating heart motions. Hence, we evaluate the performance using multiple scenarios such as tracking circular paths and AHV motions.

3.2.1. Circular path

The first set of experiments is performed in order to evaluate steering of the RADS along circular paths using the MPC strategy. Note, that no heart valve motion modelling is considered ( $r_{AHV} = 0$ ) while the RADS moves along the circular paths described by the reference signal ( $r_{d}$ ). The circle has a radius of 6 mm, while we use an articulating tip velocity of 2 mm/s. First, we evaluate a circular path without limitations on the RADS tip position. Subsequently, we limit the tip position to a maximum of 4 mm along the x-axis by integrating a state constraint to the MPC strategy (similar to (40)). By restricting the instrument motion to an allowable region, we demonstrate the ability of the MPC strategy to avoid sensitive tissue that could be present in surgery. These limitations can be integrated in the MPC strategy, allowing the controller to anticipate for the forbidden regions in order to avoid damage to sensitive tissue. In Figure 9(b), the allowed area of the RADS tip is depicted green, while the forbidden region is red.

Fig. 9.

Representative experimental model-predictive control results of the articulating tip of the robotically actuated delivery sheath during tracking of (a) circular reference path ( $r_{d}$ ) and ( $r_{AHV} =$ 0), (b) constrained circular reference path ( $r_{d}$ ) and ( $r_{AHV} =$ 0) and (c) and (d) aortic heart valve (AHV) motion trajectories according to $r_{AHV} = [f_{x} f_{y} 0]^{T}$ and $r_{d} =$ 0. The red line trajectory represents the reference path ( r) described in (25), while the blue line represents the actual path ( y $= [p_{x} p_{y} 0]^{T}$ ) followed by the articulating tip.

3.2.2. Aortic heart valve motion tracking

The second set of experiments is used to evaluate the novel designed MPC strategy to steer the RADS with a priori knowledge on the AHV motion. Although, other medical applications such as mitral valve surgery could potentially benefit from the integrated system, we demonstrate the capabilities of the MPC strategy by tracking AHV motions in the annular plane. The beating heart motions have been obtained by analyzing the pre-operative 2D ultrasound images as depicted in Figure 6. We evaluate the performance of the MPC with a priori knowledge of the AHV motion described by tracking reference ( $r_{AHV} = [f_{x} f_{y} 0]^{T}$ ) obtained by (31) and (32) and $r_{d} = 0$ . In order to preserve the dominant features required for ultrasound image segmentation, we limit the tip velocity to 20 mm/s in x- and y-axes of frame ( $Ψ_{0}$ ) by integrating state constraints according to (24). A RADS tip velocity that exceeds 20 mm/s, will cause dominant features to disappear in ultrasound images. This results in unreliable segmentation of the RADS tip position. Further, we limit the maximum tip deflection to 10 mm in x- and y-directions (frame ( $Ψ_{0}$ )) by (23).

3.3. Experimental results

The results of the experimental scenarios are reported in Table 2, while a single representative of each experiment can be found in Figure 9. In order to evaluate the tip tracking accuracy of the RADS, the experiments were repeated 12 times for circular trajectories and 30 times for AHV motions. The mean absolute errors (MAE) in the tracked tip position ( $ε_{x}$ and $ε_{y}$ ) during trajectory tracking are reported. On average, the RANSAC algorithm of the segmentation strategy completes $819 \pm 101$ iterations (single CPU core implementation). In experiments, we observe a constant delay of approximately 200 ms in obtaining the 2D ultrasound images used for feedback. Further, a compensation algorithm is used to reduce the mechanical hysteresis in the system in order to improve tip positioning accuracy.

Table 2.

Experimental results of the robotically actuated delivery sheath tip tracking for circular paths and aortic heart valve motions using the model-predictive control. The mean absolute distance error ( $ε$ ) and position errors ( $ε_{x}$ and $ε_{y}$ ) along the x- and y-axes are provided. Further, the standard deviation for $N_{r}$ repetitions is reported.

Case	$N_{r}$	$ε_{x}$ [mm]	$ε_{y}$ [mm]	$ε$ [mm]
Case	$N_{r}$	Circle	12	0.73±0.03	0.50±0.10	0.89±0.08
Constrained circle	12	0.74±0.11	0.61±0.17	0.97±0.17
Aortic heart valve	30	1.06±0.43	1.29±0.38	1.68±0.53

Experiments showed tracking of a circular trajectory with MAE of 0.73 and 0.50 mm in the x- and y-axes, respectively. The observed mean absolute distance error is 0.89 mm. A slightly higher error is observed in the x-axis compared to the y-axis. This could be attributed to the difference in mechanical properties such as friction and backlash of the two uncoupled pulley systems (each control a single DOF). Nonetheless, the results show a significant improvement in tracking performance compared to the mean positioning errors of approximately 2 mm along the x- and y-axes demonstrated in previous work (Vrooijink et al., 2014). Further, in the experiments that demonstrate the circular trajectory in which the instrument motion is limited, we observe a MAE of 0.74 and 0.61 mm in the x- and y-axes, respectively. The corresponding mean absolute distance error is 0.97 mm. The observed error correspond to the results of a circular trajectory without instrument limitations. Hence, we demonstrate that the MPC strategy is capable of avoiding areas that could potentially be sensitive tissue without degrading performance. By tracking AHV motions, we observe a MAE of 1.06 and 1.29 mm in the x- and y-axes, respectively. The observed mean absolute distance error is 1.68 mm. A higher error is observed in tracking AHV motions compared with circular trajectories. This could be explained by the fast moving AHV motions that impedes the accuracy of tracking. Further, we observe a decrease in tracking accuracy along the y-axis for AHV motions. The decrease in tracking accuracy could be explained by the instrument velocity which is limited to 20 mm/s in MPC in order to enable RADS detection in ultrasound images. The MPC strategy considers both the tip velocity limitation and AHV motion in order to optimize tracking accuracy.

Further, the ultrasound measurement delay introduced in the feedback signal of the MPC strategy is among the main contributor to the tracking error. This measurement delay is considered to be the combined result of pre-processing within the ultrasound system, transferring images to the computer using a capturing device and instrument segmentation. However, it cannot be ruled out that tracking accuracy is affected by (non-)linear friction, tendon elongation or crosstalk between the two tendon pairs (i.e. controlling one tendon pair influences the other tendon pair). Note, that we evaluate our system using pre-operative patient data which leads to simplifications with respect to heart-rate variability. However, these simplifications could potentially be eliminated by administering medicine that reduce the heart rate and apply over-pacing to effectively control the patient’s heart rate in a predictable manner. This improves the ability to anticipate the AHV motion. Nonetheless, our method demonstrates tracking of AHV motions based on pre-operative ultrasound patient data. This indicates that our system, with further development, could provide cardiac motion compensation to a wide class of cardiovascular applications performed without a heart–lung machine.

4. Conclusions and future work

In this study, we present and evaluate a novel RADS capable of autonomously and accurately compensating for AHV motions by using an MPC strategy. We develop and incorporate kinematic models of the RADS within the MPC strategy. In order to accurately evaluate the RADS tip position, we use a clinically available 2D ultrasound transducer which is orientated perpendicular to the tip of the RADS. An online segmentation algorithm is developed in order to provide feedback of the RADS tip position for the MPC strategy. Pre-operative ultrasound images of a patient are used to evaluate the AHV motion. Further, mechanical hysteresis is addressed by a compensation algorithm in order to improve the tip positioning accuracy. The novel integrated system is capable of controlling the articulating tip of the RADS in order to compensate for AHV motions. In experiments, we demonstrate evidence that the RADS tracks the AHV motions with a mean absolute distance error for AHV motions of 1.68 mm. Hence, we potentially improve and enable new applications in cardiovascular surgery performed without a heart–lung machine.

In future work, we intend to address accuracy problems introduced by measurement delay and instrument friction. Electro-magnetic instrument tip tracking will be integrated in order to reduce feedback delay and to improve robustness in control. The prototype device used in this study will be replaced by a flexible steerable catheter in order to enable applications in ablation, aortic and mitral valve surgery. Further, we plan to combine instrument and intra-operative AHV motion segmentation of 2D and 3D ultrasound images in order to improve the tracking performance in a clinically relevant scenario. We continuously aim to improve the robustness and accuracy of the integrated system. Nonetheless, the presented framework for modelling, imaging and control is applicable to a range of continuum-style robots and catheters. Our current study evaluates compensation of AHV motions using an MPC strategy. Hence, we have demonstrated the feasibility and potential for steerable catheters to compensate for cardiac motions in cardiovascular interventions such as aortic and mitral valve surgery.

Footnotes

Appendix A: Index to Multimedia Extension

Archives of IJRR multimedia extensions published prior to 2014 can be found at http://www.ijrr.org, after 2014 all videos are available on the IJRR YouTube channel at http://www.youtube.com/user/ijrrmultimedia

SAGE-Journals-Accessible-Video-Player 10.1177/0278364917691113.M1 sj-vid-1-ijr-10.1177_0278364917691113 Extension 1

This video demonstrates some representative results of tracking circular paths and aortic heart valve motions using model-predictive control.

Appendix B: Nomenclature

k	Discrete time variable	[-]
$ψ_{*}$	Pulley angle, where $* \in {x, y}$	[rad]
$δ_{*}$	Tendon displacement, where $* \in {x, y}$	[m]
$κ$	Arc curvature	[ $\frac{1}{m}$ ]
r	Arc radius	[m]
$ϕ$	Arc plane angle	[rad]
$ℓ$	Arc length	[m]
$θ$	Arc bend angle	[rad]
$l_{t}$	Rigid link length	[m]
$r_{p}$	Pulley radius	[m]
$t_{i}$	Tension wires, where $i = 1, \dots, 4$	[–]
$l_{i}$	Tendon lengths, where $i = 1, \dots, 4$	[m]
$ϕ_{i}$	Individual tendon angle, where $i = 1, \dots, 4$	[m]
$d_{b}$	Distance between backbone and tendons	[m]
$Ψ_{*}$	Coordinate system, where $* \in {0, i, t, u}$	[–]
$H_{i}^{0}$	Homogeneous transformation matrix expressing the intermediate frame in the reference frame ( $H_{i}^{0} \in SE 3$ )	[–]
$H_{t}^{i}$	Homogeneous transformation matrix expressing the tip frame in the intermediate frame ( $H_{t}^{i} \in SE 3$ )	[–]
$H_{t}^{0}$	Homogeneous transformation matrix expressing the tip frame in the reference frame ( $H_{t}^{0} \in SE 3$ )	[–]
$L_{t}^{i}$	Rigid link section by a translation along the z-axis of the intermediate frame ( $L_{t}^{i} \in R^{3}$ )	[–]
$r$ , $r_{t}^{0}$	Referenced tip position expressed in the reference frame ( $r$ , $r_{t}^{0} \in R^{3}$ )	[–]
$r_{*}$	Referenced tip position, where $* \in {x, y, z}$	[m]
$o_{t}$	Origin of the articulating tip frame ( $o_{t} \in R^{3}$ )	[–]
$p$ , $p_{t}^{0}$	Tip position expressed in the reference frame ( $p$ , $p_{t}^{0} \in R^{3}$ )	[–]
$p_{*}$	Tip position, where $* \in {x, y, z}$	[m]
$u$	Arc parameters, $u \in R^{2}$	[–]
G	Rate of convergence feedback gain	[–]
e	Estimation error	[m]
$ε_{th}$	Error tolerance	[m]
$A_{c}$	Set of detected edge points	[–]
$x_{c}$	Edge point obtained from Canny edge detector ( $x_{c} \in R^{2}$ )	[–]
$H_{3}$	Set of three randomly selected edge points	[–]
$m_{c}$	Algebraic circle model parameters	[–]
n	Number of random sampling consensus iterations	[–]
$J_{c}$	Random sampling consensus cost	[–]
$S_{c}$	Random sampling consensus set (inliers)	[–]
$R^{2}$	Data fitting discrepancy	[–]
$e$	Zero mean white noise ( $e, \hat{e} \in R^{2}$ )	[–]
$v$	Control input signal ( $v \in R^{2}$ )	[–]
y	MPC tip position ( $y \in R^{3}$ )	[–]
$p_{\min}, p_{\max}$	Inequality constraint lower and upper position bound ( $p_{\min}, p_{\max} \in R^{3}$ )	[m]
${\overset{\cdot}{p}}_{\min}, {\overset{\cdot}{p}}_{\max}$	Inequality constraint lower and upper velocity bound ( ${\overset{\cdot}{p}}_{\min}, {\overset{\cdot}{p}}_{\max} \in R^{3}$ )	[ $m / s$ ]
$T_{s}$	Sampling time	[s]
$a$	Variations in arc parameters ( $a \in R^{2}$ )	[–]
$r_{AHV}$	Aortic heart valve annulus motion reference ( $r_{AHV} \in R^{3}$ )	[–]
$r_{d}$	Desired tip position ( $r_{d} \in R^{3}$ )	[–]
J	Generalized predictive control cost	[–]
N	Prediction horizon	[–]
$N_{m}$	Minimum cost horizon	[–]
$N_{c}$	Control horizon	[–]
$λ$	Control input weighting matrix ( $λ \in R^{2 \times 2}$ )	[–]
$z$	Cost signal ( $z \in R^{3}$ )	[–]
$Γ$	Diagonal selection matrix ( $Γ \in R^{5 \times 5}$ )	[–]
$f_{*}$	AHV reference position, where $* \in {x, y}$	[m]
$a_{}$ , $b_{}$	The Fourier series coefficients	[–]
$ω$ , $ω_{*}$	The Fourier series frequencies, where $* \in {x, y}$	[m]
$S_{r}$	Sample rate Fourier series	$[\frac{1}{s}]$
$H_{r}$	Heart rate	[bpm]
$e_{rr}$	The tip position error ( $e_{rr} \in R^{3}$ )	[–]
$x$	The MPC state variable ( $x \in R^{5}$ )	[–]
$w$	The MPC external signal ( $w \in R^{6}$ )	[–]
$d$	The MPC disturbance ( $w \in R^{3}$ )	[–]
$A, B, C, D$	Multiple-input and multiple-output state space realization of the CSPCP	[–]
$υ$	Equality constraint signal ( $υ \in R^{2}$ )	[–]
$ϒ$	Equality constraint ( $ϒ \in R^{2}$ )	[–]
$χ, χ_{1}, χ_{2}$	Inequality constraint signals ( $χ_{1}, χ_{2} \in R^{5}, χ \in R^{10}$ )	[–]
$X, X_{1}, X_{2}$	Inequality constraints ( $X_{1}, X_{2} \in R^{5}, X \in R^{10}$ )	[–]
$Δ c$	Compensated angular pulley velocity ( $Δ c \in R^{2}$ )	[–]
$Δ ψ_{*}$	Angular pulley velocity, where $* \in {x, y}$	[rad/s]
$α_{*}^{+}$	Positive contact positions, where $* \in {x, y}$	[m]
$α_{*}^{-}$	Negative contact positions, where $* \in {x, y}$	[m]
$ε_{*}$	Mean absolute error in the tracked tip position, where $* \in {x, y}$	[m]
$ε$	Mean absolute tip distance error	[m]
$N_{r}$	Number of experimental repetitions	[–]

Appendix C: Concept of prediction

Based on the model description (34), the concept of prediction can be introduced according to

(51)

{\tilde{s}}_{p, (k)} = {\tilde{C}}_{p} x_{(k)} + {\tilde{D}}_{p 1} e_{(k)} + {\tilde{D}}_{p 2} {\tilde{w}}_{(k)} + {\tilde{D}}_{p 3} {\tilde{v}}_{(k)}

where any signal $s_{p, (k)}$ as presented in (34) can be used to formulate predictions at each time instant (k), that considers all signals over horizon N. The prediction of signal $s_{p, (k)}$ is given by signal vector ${\tilde{s}}_{p, (k)}$ according to

(52)

{\tilde{s}}_{p, (k)} = [\begin{matrix} {\hat{s}}_{p, (k | k)} \\ {\hat{s}}_{p, (k + 1 | k)} \\ ⋮ \\ {\hat{s}}_{p, (k + N - 1 | k)} \end{matrix}]

The prediction presented in (51) is based on future control and reference signals described by signal vectors

(53)

{\tilde{v}}_{(k)} = [\begin{matrix} {\hat{v}}_{(k | k)} \\ {\hat{v}}_{(k + 1 | k)} \\ ⋮ \\ {\hat{v}}_{(k + N - 1 | k)} \end{matrix}] and {\tilde{w}}_{(k)} = [\begin{matrix} {\hat{w}}_{(k | k)} \\ {\hat{w}}_{(k + 1 | k)} \\ ⋮ \\ {\hat{w}}_{(k + N - 1 | k)} \end{matrix}]

respectively. Note that ${\tilde{w}}_{(k)}$ contains a priori knowledge of the AHV position. We can obtain matrices ${\tilde{C}}_{p}$ , $\tilde{D_{p 1}}$ , $\tilde{D_{p 2}}$ and $\tilde{D_{p 3}}$ by using successive substitution of the state equation (33) into (34) which yields

(54)

{\tilde{C}}_{p} = [\begin{matrix} C_{p} \\ C_{p} A \\ C_{p} A^{2} \\ ⋮ \\ C_{p} A^{N - 1} \end{matrix}], {\tilde{D}}_{p 1} = [\begin{matrix} D_{p 1} \\ C_{p} B_{1} \\ C_{p} A B_{1} \\ ⋮ \\ C_{p} A^{N - 2} B_{1} \end{matrix}]

and

(55)

{\tilde{D}}_{pq} = [\begin{matrix} D_{pq} & 0 & \dots & 0 & 0 \\ C_{p} B_{q} & D_{pq} & ⋱ & ⋮ & ⋮ \\ C_{p} A B_{q} & C_{p} B_{q} & ⋱ & 0 & 0 \\ ⋮ & ⋱ & D_{pq} & 0 \\ C_{p} A^{N - 2} B_{q} & \dots & C_{p} B_{q} & D_{pq} \end{matrix}]

where subscript $q = 2, 3$ is used to indicate the corresponding matrix to complete the prediction. The concept of prediction is used to provide a prediction model of the CSPCP signals described in (34), (35) and (36).

Acknowledgements

The authors would like to thank Dr J Scheltes (DEAM Corporation, the Netherlands), Professor P Breedveld and R van Starkenburg from the Delft University of Technology, the Netherlands.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research is supported by the Dutch Technology Foundation STW (iMIT-Instruments for Minimally Invasive Techniques Interactive: Multi-Interventional Tools (Project: MULTI)), which is part of the Netherlands Organisation for Scientific Research (NWO) and partly funded by the Ministry of Economic Affairs, Agriculture and Innovation.

References

Aldrich

(2007) Basic physics of ultrasound imaging. Critical Care Medicine 35(5): S131–S137.

Alfieri

Bonis

Maisano

Canna

(2007) Future directions in degenerative mitral valve repair. Seminars in Thoracic and Cardiovascular Surgery 19(2): 127–132.

Alwan

(2011) Global status report on noncommunicable diseases 2010. Geneva: World Health Organization.

Bar-Shalom

Kirubarajan

(2001) Estimation with Applications to Tracking and Navigation. New York: John Wiley & Sons, Inc.

Bardou

Zanne

Nageotte

De Mathelin

(2010) Control of a multiple sections flexible endoscopic system. In: Proceedings of the IEEE International Conference on Intelligent Robots and Systems (IROS), Taipei, Taiwan, pp. 2345–2350.

Bebek

Cavusoglu

(2007) Intelligent control algorithms for robotic-assisted beating heart surgery. IEEE Transactions on Robotics 23(3): 468–480.

Bowthorpe

Alvarez Garcia

Tavakoli

(2014) GPC-based teleoperation for delay compensation and disturbance rejection in image-guided beating-heart surgery. In: Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), Hong Kong, China, pp. 4875–4880.

Breedveld

(2010) Steerable laparoscopic cable-ring forceps. In: Proceedings of the Design of Medical Devices Conference (DMD), Minneapolis, USA, p. 7.

Camarillo

Milne

Carlson

Zinn

Salisbury

(2008) Mechanics modeling of tendon-driven continuum manipulators. IEEE Transactions on Robotics 24(6): 1262–1273.

10.

Cribier

Eltchaninoff

Tron

. (2004) Early experience with percutaneous transcatheter implantation of heart valve prosthesis for the treatment of end-stage inoperable patients with calcific aortic stenosis. Journal of the American College of Cardiology 43(4): 698–703.

11.

Ding

Goldman

Allen

Fowler

Simaan

(2010) Design, simulation and evaluation of kinematic alternatives for insertable robotic effectors platforms in single port access surgery. In: Proceedings of the IEEE International Conference on Robotics and Automation (ICRA). Anchorage, USA, pp. 1053–1058. DOI:10.1109/ROBOT.2010.5509244.

12.

Dupont

Lock

Itkowitz

Butler

(2010) Design and control of concentric-tube robots. IEEE Transactions on Robotics 26(2): 209–225.

13.

Falk

(2002) Manual control and tracking - a human factor analysis relevant for beating heart surgery. The Annals of Thoracic Surgery 74(2): 624–628.

14.

Ferreau

Kirches

Potschka

Bock

Diehl

(2014) qpOASES: A parametric active-set algorithm for quadratic programming. Mathematical Programming Computation 6(4): 327–363.

15.

Forsyth

Ponce

(2002) Computer vision a modern approach. Englewood Cliffs, NJ: Prentice Hall Professional Technical Reference.

16.

Gangloff

Ginhoux

de Mathelin

Soler

Marescaux

(2006) Model predictive control for compensation of cyclic organ motions in teleoperated laparoscopic surgery. IEEE Transactions on Control Systems Technology 14(2): 235–246.

17.

Ginhoux

Gangloff

de Mathelin

Soler

Sanchez

Marescaux

(2005) Active filtering of physiological motion in robotized surgery using predictive control. IEEE Transactions on Robotics 21(1): 67–79.

18.

Gosline

Vasilyev

Butler

. (2012) Percutaneous intracardiac beating-heart surgery using metal MEMS tissue approximation tools. The International Journal of Robotics Research 31(9): 1081-1093.

19.

Himbert

Descoutures

Al-Attar

. (2009) Results of transfemoral or transapical aortic valve implantation following a uniform assessment in high-risk patients with aortic stenosis. Journal of the American College of Cardiology 54(4): 303–311.

20.

Hong

Dohi

Hashizume

Konishi

Hata

(2004) An ultrasound-driven needle-insertion robot for percutaneous cholecystostomy. Physics in Medicine and Biology 49(3): 441–455.

21.

Jayender

Azizian

Patel

(2008) Autonomous image-guided robot-assisted active catheter insertion. IEEE Transactions on Robotics 24(4): 858–871.

22.

Johansson

Nozohoor

Kimblad

Harnek

Olivecrona

Sjögren

(2011) Transapical versus transfemoral aortic valve implantation: A comparison of survival and safety. The Annals of Thoracic Surgery 91(1): 57–63.

23.

Kesner

Howe

(2014) Robotic catheter cardiac ablation combining ultrasound guidance and force control. The International Journal of Robotics Research 33(4): 631–644.

24.

Kinnaert

(1989) Adaptive generalized predictive controller for MIMO systems. International Journal of Control 50(1): 161–172.

25.

Koolwal

Barbagli

Carlson

Liang

(2010) An ultrasound-based localization algorithm for catheter ablation guidance in the left atrium. The International Journal of Robotics Research 29(6): 643–665.

26.

Nadeau

Ren

Krupa

Dupont

(2015) Intensity-based visual servoing for instrument and tissue tracking in 3D ultrasound volumes. IEEE Transactions on Automation Science and Engineering 12(1): 367–371.

27.

Nakamura

Kishi

Kawakami

(2001) Heartbeat synchronization for robotic cardiac surgery. In: Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), volume 2, Seoul, Korea, pp. 2014–2019.

28.

Neshat

HRS

Patel

(2008) Real-time parametric curved needle segmentation in 3D ultrasound images. In: Proceedings of the IEEE RAS EMBS International Conference on Biomedical Robotics and Biomechatronics (BioRob), Scottsdale, USA, pp. 670–675.

29.

Neubach

Shoham

(2010) Ultrasound-guided robot for flexible needle steering. IEEE Transactions on Biomedical Engineering 57(4): 799–805.

30.

Novotny

Stoll

Vasilyev

. (2007) GPU based real-time instrument tracking with three-dimensional ultrasound. Medical Image Analysis 11(5): 458–464.

31.

Ortmaier

Groger

Boehm

Falk

Hirzinger

(2005) Motion estimation in beating heart surgery. IEEE Transactions on Biomedical Engineering 52(10): 1729–1740.

32.

Reilink

Stramigioli

Misra

(2013) Image-based hysteresis reduction for the control of flexible endoscopic instruments. Mechatronics 23(6): 652–658.

33.

Reiter

Bajo

Iliopoulos

Simaan

Allen

(2012) Learning-based configuration estimation of a multi-segment continuum robot. In: Proceedings of the IEEE International Conference on Biomedical Robotics and Biomechatronics (BioRob), Rome, Italy, pp. 829–834.

34.

Richa

Bó

Poignet

(2010) Beating heart motion prediction for robust visual tracking. In: Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), Anchorage, USA, pp. 4579–4584.

35.

Rone

Ben-Tzvi

(2014) Continuum robot dynamics utilizing the principle of virtual power. IEEE Transactions on Robotics 30(1): 275–287.

36.

Seeburger

Borger

Tschernich

. (2010) Transapical beating heart mitral valve repair. Circulation: Cardiovascular Interventions 3(6): 611–612.

37.

Stoll

Ren

Dupont

(2012) Passive markers for tracking surgical instruments in real-time 3D ultrasound imaging. IEEE Transactions on Medical Imaging 31(3): 563–575.

38.

Tamburino

Ussia

Maisano

. (2010) Percutaneous mitral valve repair with the mitraclip system: acute results from a real world setting. European Heart Journal 31(11): 1382–1389.

39.

Tuna

Franke

Bebek

Shiose

Fukamachi

Cavusoglu

(2013) Heart motion prediction based on adaptive estimation algorithms for robotic-assisted beating heart surgery. IEEE Transactions on Robotics 29(1): 261–276.

40.

van den Boom

Stoorvogel

(2012) Model predictive control. Lecture Notes for the Dutch Institute of Systems and Control (DISC), The Netherlands.

41.

Vrooijink

Abayazid

Misra

(2013) Real-time three-dimensional flexible needle tracking using two-dimensional ultrasound. In: Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), Karlsruhe, Germany, pp. 1680–1685.

42.

Vrooijink

Ellenbroek

TTM

Breedveld

Grandjean

Misra

(2014) A preliminary study on using a robotically-actuated delivery sheath (RADS) for transapical aortic valve implantation. In: Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), Hong Kong, China, pp. 4380–4386.

43.

Walther

Dewey

Borger

. (2009) Transapical aortic valve implantation: Step by step. The Annals of Thoracic Surgery 87(1): 276–283.

44.

Webster

III Jones

(2010) Design and kinematic modeling of constant curvature continuum robots: A review. The International Journal of Robotics Research 29(13): 1661–1683.

45.

Wong

Cheung

Webb

Carere

Lichtenstein

(2010) Technical considerations to avoid pitfalls during transapical aortic valve implantation. The Journal of Thoracic and Cardiovascular Surgery 140(1): 196–202.

46.

Cheung

Lichtenstein

. (2010) Transapical transcatheter aortic valve implantation: Follow-up to 3 years. The Journal of Thoracic and Cardiovascular Surgery 139(5): 1107–1113.

47.

Yuen

Novotny

Howe

(2008) Quasiperiodic predictive filtering for robot-assisted beating heart surgery. In: Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), Pasadena, USA, pp. 3875–3880.

Model predictive control of a robotically actuated delivery sheath for beating heart compensation

Abstract

Keywords

1. Introduction

1.1. Related work

1.2. Contributions

2. Methods

2.1. Instrument modelling

2.1.1. Forward kinematics

2.1.2. Inverse kinematics

2.1.3. Differential kinematics

2.2. Ultrasound image segmentation

2.3. Model predictive control

2.3.1. Model description

2.3.2. Control objective

2.3.3. AHV modelling

2.3.4. Controller design

2.3.5. Hysteresis compensation

3. Experiments

3.1. Experimental setup and materials

3.2. Experimental scenarios

3.2.1. Circular path

3.2.2. Aortic heart valve motion tracking

3.3. Experimental results

4. Conclusions and future work

Footnotes

Appendix A: Index to Multimedia Extension

Appendix B: Nomenclature

Appendix C: Concept of prediction

Acknowledgements

Funding

References