Evaluation of an Active Disturbance Rejection Controller for Ophthalmic Robots with Piezo-driven Injector

1. Introduction

Retinal vein occlusion (RVO) occurs when a blockage of one or more small veins carries blood away from the retina. RVO has been reported to be the second most common retinal vascular occlusive disease[1]. Its impact on vision can be severe, causing blurred and distorted vision[2]. Currently, there is no available treatment to resolve RVO directly. A possible treatment is retinal vein cannulation (RVC), which involves the insertion of a cannulation needle and light pipe through the scleral ports (diameter: about 0.6mm) to inject a clot-dissolving drug into the retinal vein. Figure 1 gives a visual representation of the procedure of RVC.

Performing micro-surgical tasks in RVC requires highly specialized skills that are at the limits of human motor capabilities. The surgeons’ hand tremor amplitudes exceeds 150 μm[3], which is significantly larger than the diameters of the retinal veins (50-150 μm[4]). The surgeons are not able to perceive the contact force at the needle tip, which is below 10mN[5]. The limited motor and force control, degraded manual dexterity, and possible visual, physical, and mental fatigue of a surgeon may intensify their hand tremors and increase the risk of sight-damaging iatrogenic trauma to the eye.

A technically feasible remedy to these challenges is robot-assisted RVC. Figure 1 illustrates the three processes involved in robot-assisted RVC: (1) Approach – where the robot drives the needle approaches the retinal vein. (2) Contact – where the robot ensures that the contact force reaches the expected value, making a puncture into the vein without causing any damage to the retina. (3) Inject- where the robot keeps the needle steady and injects the clot-dissolving drug into the retina vein.

Current robots designed for robotic assisted vascular catheterization (RVC) have primarily focused on improving the precision of surgical tools and filtering out the surgeon’s hand tremors. For instance, leveraging the natural steady motion of robotic arms [6,7,8], compensating the tool tip’s motion with hand-held devices [9], wireless magnetic control of an untethered microrobot [10], localization of the tool tip with OCT[11]/computer vision [12], sensing tool-tissue proximity with bio-impedance [13]. While precise needle movement is the foundation of robot-assisted RVC, the lack of a contact force controller presents a safety challenge. A robotic system which is capable of sensing and controlling the contact force is advantageous.

Researchers from Johns Hopkins University have developed force-sensing tools that utilize fiber Bragg grating (FBG) sensors. With these sensors, the tools can measure the contact force, scleral force, and insertion depth. The force measurement can be integrated into the robot control law [14,15] or provide auditory substitution [16]. Dual force constraint controller[17], hybrid position/force controller [18] are also proposed to keep the force within desired range. In the contact process, the robot is expected to track the desired force during the contact process, not just limit the contact force.

Ebrahimi et al [19] model the sclera with linear elastic contact model and propose an adaptive controller to track the sclera force. However, the stiffness fail to describe the retina vein’s viscoelastic characteristics. Additionally, the viscoelastic parameters of eye tissue are typically related to a variety of factors, such as age, gender, and species (human or different experimental animals) [20], making the contact force highly uncertain and complex. To improve the force control performance, a viscoelastic contact model and a robust controller are necessary.

The parameter uncertainty could be regarded as inner disturbance in viscoelastic contact model. Active disturbance rejection control (ADRC), presented by Han[21], is a control strategy with strong adaptability, robustness, and independence on the accurate parameters of the system model. ADRC is proved to have good rejection ability against the disturbance and noise[22,23]. To the best of our knowledge, the ADRC has not been applied to reject the inner disturbance in viscoelastic contact model and control the “needle-retina vein” contact force.

This preliminary study establishes a contact model based on viscoelastic theory to capture the nonlinear, uncertain behavior of vessel. Then, on the basis of the viscoelastic contact model, an active disturbance rejection force controller (ADRC) is designed to reject the disturbance caused by the parameter uncertainty, and precisely control the contact force. Section 2 introduces the proposed ophthalmic surgical robot and builds the viscoelastic contact model. Section 3 focuse on the design of the ADRC. Section 4 analyzes the disturbance rejection ability through simulation, while Section 5 presents the experimental setup and results. Section 5 concludes this paper.

2. Materials and Methods

2.1. Ophthalmic Surgical Robot

Figure 2a and Figure 2b give the components and mechanism configuration of the ophthalmic surgical robot system, respectively. The proposed robot can perform the instruments’ orientation, insertion, and rotation. The instruments’ orientation are adjusted by four linear motors (LM1-4, MSR080, Paker, USA). We can rotate the instrument around or along axis y_B by changing the stroke of LM1 and LM3. Similarly, we can rotate the instrument around or along axis x_B by changing the stroke of LM2 and LM4. The proposed ophthalmic surgical robot can also perform continuous curvilinear capsulorhexis or membrane peeling with forceps [8,25].

The piezo-driven injector consists of a piezoelectric linear stage (CONEX-SAG-LS48P, NewPort, USA, P₅), a rotation motor (DCX10S, Maxson, Switzerland, R₅) and the instrument shown in Figure 2c. The positioning accuracy of the linear motors and piezoelectric linear stage are 0.01 mm and 25 nm, respectively. The detailed components of instrument is given in Figure 2c. The instrument consists of a force sensor, a connector, a Luer fitting and a cannulation needle (41G). The needle is also equipped with three FBG sensors. The force sensor (LSB200-20g, Futek, USA, precision: 0.2 mN) aims to measure the axial force (z, Figure 2b-c) acting at the needle.

2.2. Viscoelastic Contact Model

Figure 3 illustrates the viscoelastic contact model, which models the contact force f(t) with two parallel Maxwell viscoelastic models [26]. Each Maxwell model combines a spring with a damper in serial.

The mathematical description of the viscoelastic contact model is given in Eq.(1).

$$f\left(t\right)=\left(E_1+E_2\right)\frac{dx\left(t\right)}{dt}-\left({\tau }_1+{\tau }_2\right)\frac{df\left(t\right)}{dt}$$

(1)

where x(t) is the axial movement of PLS1, τ_i = η_i/E_i, i=1,2. E and η represent the spring stiffness constant and the damping factor, respectively.

E and τ are the mechanical parameters of the retinal vein and could be written as below:

$E=E_{ave}+\Delta E,\tau ={\tau }_{ave}+\Delta \tau ,$ E=E_ave+△E, τ=τ_ave+△τ

where ave means the mean value measured by experiments, △E and △τ are the uncertainty caused by many factors, such as age, gender, experiment models (in-vivo/ex-vivo, pig/rabbit/human),etc.

Eq.1 could be rewritten in frequency domain:

$$F\left(s\right)=\left(\frac{E_1{\tau }_{1}s}{1+{\tau }_{1}s}+\frac{E_2{\tau }_{2}s}{1+{\tau }_{2}s}\right)X\left(s\right)$$

(2)

where s is the Laplace operator.

The step response of Eq.2 in time-domain is shown as below:

$$f\left(t\right)=E_1{e}^{-t/{\tau }_1}+E_2{e}^{-t/{\tau }_2}$$

(3)

Eq.3 is used to identify E_i,ave, τ_i,ave, i=1,2 in contact model, which will be given in Section 5.

△E and △η are the main disturbance of contact model, which makes it hard to control f(t) precisely. In this paper, a LADRC-based force controller with disturbance-rejection ability is first presented and designed for viscoelastic contact model.

2.3. Active Disturbance Rejection Force Controller

Figure 4 gives the block diagram of the ADRC. The ADRC consists of 3 parts: linear extended state observer (LESO), linear state error feedback (LSEF), and a numerical calculation module based on the fourth-order Runge Kutta (RK4) method. This part gives the detail design of the controller.

First, we rewrite Eq.2 into Eq.4

$$\frac{F\left(s\right)}{X\left(s\right)}=\frac{b_{1}s}{1+a_{1}s}+\frac{b_{2}s}{1+a_{2}s}$$

(4)

where, ${\eta }_i=E_i{\tau }_i$, $a_i={\tau }_i,b_i={\eta }_i,$i=1,2

Then, we can get Eq.5 by dividing s on both sides of the Eq.4.

$$\frac{Y\left(s\right)}{U\left(s\right)}=\frac{F\left(s\right)}{v\left(s\right)}=\frac{b_1}{1+a_{1}s}+\frac{b_2}{1+a_{2}s}$$

(5)

where Y(s) is regarded as the output of viscoelastic contact model. v(s) is the velocity of the PLS1, and it is defined as the input U(s).

Eq.5 could be rewritten as below:

$$\frac{Y\left(s\right)}{U\left(s\right)}=\frac{\left(b_1+b_2\right)+\left(a_1{b}_2+a_2{b}_1\right)s}{1+\left(a_1+a_2\right)s+a_1{a}_2{s}^2}=\frac{B_{1}s+B_2}{s^2+A_{2}s+A_3}$$

(6)

Next, we define the median parameter Y₁(s) in Eq.7

$$Y_1\left(s\right)=\frac{1}{s^2+A_{2}s+A_3}U\left(s\right)$$

(7)

Submitting Eq.7 into Eq.6, we get Eq.8

$$Y\left(s\right)=Y_1\left(s\right)\left(B_{1}s+B_2\right)$$

(8)

We get Eq.9-10 by rewriting Eq.7-8 in time-domain. Eq.9 is the state equation, Eq.10 is the observation equation.

$${\ddot{y}}_1=u-A_2{\dot{y}}_1-A_3{y}_1$$

(9)

$$y=B_1{\dot{y}}_1+B_2{y}_1$$

The uncertainty (△E and △τ) exists in both state equation and observation equation. In this paper, we only deal with the uncertainty in state equation. Then, we obtain Eq.11 by representing the uncertainty (△E and △τ) with ω.

$$\begin{array}{l}{\ddot{y}}_1=f+b_{0}u\\ y=B_1{\dot{y}}_1+B_2{y}_1\\ f=-A_2{\dot{y}}_1-A_3{y}_1+\omega +\left(1-b_0\right)u\end{array}$$

(10)

where b₀ is the known part in u, f is the sum of uncertainty and unknown disturbance.

Next, we define the state variable $x={\left[\begin{array}{ccc}{y}_1& {\dot{y}}_1& f\end{array}\right]}^{\mathrm{T}}$ and rewritte Eq.11 into Eq.12.

$$\begin{array}{l}x=Ax+Bu+E\dot{f}\\ y=Cu\end{array}$$

(11)

where,

$$A=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right],B=\left[\begin{array}{c}0\\ b_0\\ 0\end{array}\right],E=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],C=\left[\begin{array}{ccc}{B}_2& B_1& 0\end{array}\right]$$

Then, the LESO could be obtained as below:

$$\begin{array}{l}\dot{z}=Az+Bu+L\left(y-Cz\right)\\ \widehat{y}=Cz\end{array}$$

(12)

where z is the state variable in LESO and aims to observe x.

Then, we rewritten Eq.13 into Eq.14 by ignoring the differentiation of f, since it could be regarded as a part of disturbance.

$$\begin{array}{l}\dot{z}=\left(A-LC\right)z+\left[B,L\right]u_c\\ y_c=z\end{array}$$

(13)

where, $u_c={\left[\begin{array}{cc}u& y\end{array}\right]}^{\mathrm{T}}$, y_c is the output of the LESO.

In Eq.14, L is the gain matrix of the LESO and defined as below:

$$L={\left[\begin{array}{ccc}3{\omega }_0& 3{\omega }_0^2& {\omega }_0^3\end{array}\right]}^T$$

(14)

where, ω₀ is the bandwidth of the LESO.

In this paper, we solve Eq.14 with RK4. Eq.15 gives the formula of the LSEF.

$$\begin{array}{l}{u}_0=k_p\left(r-z_1\right)-k_d{z}_2\\ u=\frac{u_0-z_3}{b_0}\end{array}$$

(15)

where $k_p={\omega }_c^2,k_d=2{\omega }_c$. z₁, z₂ and z₃ are obtained from the LESO, ω_c is the bandwidth of the LSEF, r is the desired force, u is the control value. u is the velocity of PLS1.

3. Disturbance Rejection Capability Analysis

This part analyzes the disturbance-rejection capability of the proposed controller. First, we get Eq.16 by expanding Eq.13.

$$\{\begin{array}{l}{\dot{z}}_1=z_2-3{\omega }_0\left(z_1-y\right)\hfill \\ {\dot{z}}_2=z_3-3{\omega }_0^2\left(z_1-y\right)+b_{0}u\hfill \\ {\dot{z}}_3=-{\omega }_0^3\left(z_1-y\right)\hfill \end{array}$$

(16)

Then, we obtain Eq.17 by rewriting Eq.16 in frequency domain.

$$\begin{array}{l}{z}_1=\frac{3{\omega }_0{s}^2+3{\omega }_0^{2}s+{\omega }_0^3}{{\left(s+{\omega }_0\right)}^3}y+\frac{b_{0}s}{{\left(s+{\omega }_0\right)}^3}u\\ z_2=\frac{\left(3{\omega }_0^{2}s+{\omega }_0^3\right)s}{{\left(s+{\omega }_0\right)}^3}y+\frac{b_{0}s\left(s+3{\omega }_0\right)}{{\left(s+{\omega }_0\right)}^3}u\\ z_3=\frac{{\omega }_0^{3}s}{{\left(s+{\omega }_0\right)}^3}y-\frac{b_0{\omega }_0^3}{{\left(s+{\omega }_0\right)}^3}u\end{array}$$

(17)

Submitting the LSEF into Eq.17, we get Eq.18

$$U\left(s\right)=\frac{1}{b_0}{G}_1\left(s\right)\left[{\omega }_c^{2}R\left(s\right)-H\left(s\right)Y\left(s\right)\right]$$

(18)

where,

$\begin{array}{c}{G}_1\left(s\right)=\frac{{\left(s+{\omega }_0\right)}^3}{{\left(s+{\omega }_0\right)}^3+2{\omega }_c{s}^2+\left({\omega }_c^2+6{\omega }_c{\omega }_0\right)s-{\omega }_0^3}\\ H\left(s\right)=\frac{\left(3{\omega }_c^2{\omega }_0+6{\omega }_c{\omega }_0^2+{\omega }_0^3\right)s^2+{\omega }_c^2{\omega }_0^3}{{\left(s+{\omega }_0\right)}^3}+\frac{\left(3{\omega }_c^2{\omega }_0^2+2{\omega }_c{\omega }_0^3\right)s}{{\left(s+{\omega }_0\right)}^3}\end{array}$

$Y\left(s\right)=\frac{1}{s^2}\left(B_2+B_{1}s\right)\left(f+b_{0}U\right)$

Then, we get Eq.19 by rewriting Eq.18.

$$Y\left(s\right)=D_r\left(s\right)R\left(s\right)+D_f\left(s\right)f\left(s\right)$$

(19)

where, D_r(s) is the input sensitivity function. D_f(s) is the disturbance sensitivity function.

Eq.20 gives the formula of D_f(s).

$$D_f\left(s\right)=\frac{k_4{s}^4+k_3{s}^3+k_2{s}^2+k_{1}s}{m_5{s}^5+m_4{s}^4+m_3{s}^3+m_2{s}^2+m_{1}s+m_0}$$

(20)

where,

$k_4=B_1$, $k_3=B_2+5B_1{\omega }_0$, $k_2=4B_1{\omega }_0^2+5B_2{\omega }_0+6B_1{\omega }_0{\omega }_c$, $k_1=B_2{\omega }_0\left(4{\omega }_0+6{\omega }_c\right)$,

$m_5=1$, $m_4=2{\omega }_c+3{\omega }_0$, $m_0=B_2{\omega }_c^2{\omega }_0^3$, $m_3=3B_1{\omega }_c^2{\omega }_0+{\omega }_c^2+6B_1{\omega }_0^2{\omega }_c+6{\omega }_c{\omega }_0+B_1{\omega }_0^3+3{\omega }_0^2$

$m_2=3B_1{\omega }_c^2{\omega }_0^2+3B_2{\omega }_c^2{\omega }_0+2B_1{\omega }_c{\omega }_0^3+6B_2{\omega }_c{\omega }_0^2+B_2{\omega }_0^3$

$m_1=B_1{\omega }_c^2{\omega }_0^3+3B_2{\omega }_c^2{\omega }_0^2+2B_2{\omega }_c{\omega }_0^3$

Next, we simulate the tracking performance of D_f(s). All simulation are performed by MATLAB. As shown in Section 2, the disturbance comes from mechanical parameters of retinal vein. In simulation process, ω_c and ω₀ are set as 16 and 60, respectively.

The range of mechanical parameters are set as below:

$E_1\in \left[0.1,2.0\right],E_2\in \left[0.01,10\right],{\tau }_1\in \left[0.1,2.0\right],{\tau }_2\in \left[0.01,10\right]$As shown in Figure 5, the value of D_f(s)f(s) converges to zero within just one second. Changes in mechanical parameters have little effect on the simulation results. The simulation results is consistent with [22]. The proposed controller is highly resistant to disturbance.

4. Experimental Results and Discussion

4.1. Experimental Setup

Figure 6 gives the experimental setup. All experiments are performed with the proposed ophthalmic surgical robot (Section 2). A force sensor (LSB200-20g, Futek, USA, precision: 0.2 mN) and a surgical cannulation needle (41G, Incyto, South Korea, diameter: 90μm) are attached to the end of the proposed robot.

Experiments using the chicken chorioallantoic membrane of 13 days fertilized chicken embryos (Figure 7) are performed to estimate the proposed force controller. All experiments are supervised by an expert surgeon with more than 20 years clinical experience. The surgeon make sure that the cannulation needle contact the vessel precisely.

The cannulation needle is driven by the piezoelectric linear stage of the proposed robot. Data is collected using an industrial computer (ARK3523, Advantech, China), which also calculates velocity command based on ADRC and sends it to the controller of piezoelectric linear stage via USB communication. The needle's moving accuracy is approximately 25 nm. Sampling and control frequency are 200 Hz.

In this article, we perform three groups of experiments to estimate the force controller. The detail setup of each group is given as below:

I. Parameter identification. This group aims to identify the vein’s mechanical parameter (E_ave and τ_ave) from step response data. In this group, the pose of robot is chosen randomly. Five embryos are used in this group. For each embryo, we adjust the robot to contact the vein, move 0.5mm along axis z, keep steady and record the contact force.

For each embryo, we choose 1 contact point with 10 replicates each point. The contact point is carefully selected by an ophthalmic doctor. Then, the model parameters (E_i,ave, τ_i,ave, i=1,2) in viscoelastic contact model can be obtained by fitting the collected force data to Eq.3.

II. Step test. This group uses another five embryos in group I. For each embryo, we adjust the robot to contact the vein, move the needle along axis z to reach the desired force. The desired force is set as 5mN. For each embryo, the pose of robot is set as 3 conditions: θ₁=45°; θ₁=30°; θ₁=0°. For each pose, we choose 3 contact points with 5 replicates each point.

III. Precise adjustment. Another five embryos are used in this group. In this group, the pose of robot is same with group II.

For each pose, we choose 2 contact points with 5 replicates each point. In each test, we first set the desired force between 1mN and 5mN. Then, we increases / decreases the desired force by 1mN or 1.5 mN every 5 seconds. In Group II and Group III, after tunning, b₀, ω_c and ω₀ are set as 20, 16 and 60, respectively.

4.2. Results and Discussion

This part gives the experimental results. Figure 8 and Table 1 show the results of Group I and Group II. Figure 8(a)-(e) are the step response curve of every embryo. The fitted parameters (E_i,ave, τ_i,ave, i=1,2) are given in Table 1. Figure 8(f)-(j) are the results of Group II.

The steady state error (SSE) and response time are also given in Table 1. As shown in Figure 8 and Table 1, for all tested embryos, the coefficient of determination R² is over 0.97. That means all the fitted parameters are reliable. But the value of E_1,ave, τ_1,ave, E_2,ave, τ_2,ave reach variation of over 81.2%, 73.4%, 403.8%, 172.3%, respectively. The model parameters do show large interindividual differences. With the proposed controller, the SSE is [0.24, 0.35] mN. The response time ranges from 1.41s to 2.23s. No obvious overshot occurs. The proposed controller is capable of rejecting the disturbance caused by the uncertainty of contact model.

As shown in Figure 9(a)-(e), the proposed controller can adjust the contact force at 1-mN interval. No obvious overshot occurs in Group III. The root mean square (RMS) value of SSE of Group II and Group III are 0.25mN, and 0.41mN, respectively. The response time ranges from 1.04s to 2.51s, the RMS and maximum SSE are 0.32 mN and 0.47 mN, respectively.

Experimental results indicate the proposed ADRC controller is capable of rejecting the disturbance and controlling the contact force precisely. The RMS value of SSE is lower than the acceptable overshot in previous works (0.9 mN) [15,17], close to the force sensing accuracy in previous work (about 0.3 mN) [7,14,18]. In Group III, the larger error associated with “decreasing the desired force” might be attributed to the hysteresis influence of vein.

About 80% of the vein puncture force is less than 5mN[5]. Hence, the same value is adopted as the desired force in Group II. The standard deviation of vein puncture force is 1mN[5]. Therefore, the same value is set as the force intervals.

In this article, the contact force is measured by a commercial force sensor, which is too large to enter the scleral port. This issue could be relieved by replacing the force sensor with customized micro-force sensor developed by FBG, which is a part of our ongoing study. The force sensor’s precision (0.2 mN) is higher than FBG force sensors (2 mN[27]). As a result, we believe the measured force data is reliable, and the proposed controller is suitable for retina vein cannulation.

It is possible to reduce the response time by increasing ω_c,but this can result in severe overshot, which increase surgical risk. Therefore, we choose to prioritize reducing the overshot, even if it means sacrificing response time.

We did not perform dynamic tracking experiments for the following reason. The proposed controller aims to reject disturbance and control the force during the contact process shown in Figure 1. The expected force should precisely puncture the vein without damaging retina. In this context, the expected force is more likely to be static and specific, rather than dynamic. Therefore, we choose to perform step response experiment and precise adjustment experiment to evaluate the disturbance rejection ability of the ADRC.

While dynamic force might have potential benefits in RVC, this paper focuses on the disturbance rejection ability of the proposed controller. The effect of dynamic force, validating this work with more samples (in-vivo/ex-vivo, pig/rabbit) belong to our ongoing study.

In our study, we use two sets of embryos, which we refer to as Group I and Group II. The reason for this is that in experimental studies, we are able to collect data and identify the viscoelastic contact model multiple times. However, in clinical settings, the parameters of each tested sample can vary significantly from the identified model, resulting in increased control error.

Our proposed method is able to accurately identify the precise value of each tested sample based on previously collected data, allowing us to control the contact force with great precision. This precise control of the contact force is the key feature of our proposed method.

The embryos in Group I are used to collect data, while the embryos in Group II are the samples tested in the clinical setting. By using the data collected from Group I, our proposed method is able to accurately control the contact force in Group II, even though the embryos in Group II are completely different from those in Group I. For this work, the ADRC is designed such that the robot can control the contact force automatically. In the future, a hybrid control method will be investigated for injecting clot-dissolving drugs, which will simultaneously limit the needle's position and contact force. This could be a potential solution to free surgeons from extremely difficult injection procedures and improve overall surgical efficiency.

5. Conclusion

This preliminary study presents an ADRC to reject the disturbance caused by vessel’s uncertain viscoelastic characteristic, and control the “needle-vessel” contact force precisely. We first model the nonlinear, uncertain behavior of the vessel with a viscoelastic contact model. Then, an ADRC is proposed based on the viscoelastic contact model. Fertilized chicken embryos are used to identify the viscoelastic model, and estimate the ADRC. Results show that ADRC can precisely control the contact force when viscoelastic parameters changes greatly. The contact force can be adjusted at 1mN interval. The RMS and maximum value of SSE are 0.32 mN and 0.41 mN, respectively. No obvious overshot occurs during experiments. This study provides a potential way to perform technically difficult RVC procedures.