Soft-NeuroAdapt: A 3-DOF Neuro-Adaptive Patient Pose Correction System For Frameless and Maskless Cancer Radiotherapy

TL;DR

This paper presents Soft-NeuroAdapt, a neuro-adaptive soft-robot system that improves patient positioning accuracy in cancer radiotherapy by compensating for head deviations using a neural network-based controller.

Contribution

It introduces a novel 3-DOF soft-robot with an adaptive neuro-controller for real-time correction of patient head deviations during radiotherapy.

Findings

Effective correction of intrafractional deviations demonstrated

Neural network controller adapts to patient movement

Ensures patient safety and comfort during treatment

Abstract

Precise patient positioning is fundamental to successful removal of malignant tumors during treatment of head and neck cancers. Errors in patient positioning have been known to damage critical organs and cause complications. To better address issues of patient positioning and motion, we introduce a 3-DOF neuro-adaptive soft-robot, called Soft-NeuroAdapt to correct deviations along 3 axes. The robot consists of inflatable air bladders that adaptively control head deviations from target while ensuring patient safety and comfort. The adaptive-neuro controller combines a state feedback component, a feedforward regulator, and a neural network that ensures correct adaptation. States are measured by a 3D vision system. We validate Soft-NeuroAdapt on a 3D printed head-and-neck dummy, and demonstrate that the controller provides adaptive actuation that compensates for intrafractional deviations in patient positioning.

Figures

Tables

Table 1. TABLE I: Plane Segmentation Algorithm

1. for i = 1 to N do

2.   sample non-collinear points {$p_i,p_j,p_k$} from $𝒫$

3.   calculate the model coefficients $a\text{x}+b\text{y}+c\text{z}=d$

4.   find distances from all $\text{p}\in 𝒫$ to the plane $(a,b,c,d)$

5.   store points ${\text{p}}^{\ast }\in 𝒫$ that satisfy the model hypothesis,

$0\le \left|d\right|\le |d_{max}$.

6. return maximum of the stored points ${\text{p}}^{\ast }$.

Equations

$$\begin{array}{r}x\le \text{min}\Vert {\text{p}}_i-{\text{p}}_j{\Vert }_2\end{array}$$

$$\text{Q}({\mathbf{\Sigma }}_{px})=\left[\begin{array}{cc}\text{tr}({\mathbf{\Sigma }}_{px})& {\mathbf{\Delta }}^T\\ \mathbf{\Delta }& {\mathbf{\Sigma }}_{px}+{\mathbf{\Sigma }}_{px}^T-\text{tr}({\mathbf{\Sigma }}_{px}){\text{I}}_3\end{array}\right].$$

$${\overrightarrow{q}}_T={\overrightarrow{\mu }}_x-\text{R}({\overrightarrow{q}}_R){\overrightarrow{\mu }}_p$$

$$\begin{array}{rl}{Z}^N=\{& u(k),u(k-1),\dots u(k-n_u),\\ & y(k),y(k-1),\dots ,y(k-n_y)\}\end{array}$$

$$\dot{\text{y}}=\text{A}\text{y}+\text{B}\mathbf{\Lambda }\left(\text{u}-f(\text{y, u})\right)+\text{w}(k)$$

$${\dot{\text{y}}}_m={\text{A}}_m{y}_m+{\text{B}}_m\text{r},$$

$$\text{u}={\hat{\text{K}}}_y^{\text{T}}\text{y}+{\hat{\text{K}}}_r^{\text{T}}\text{r}+\hat{f}(\text{y},\text{u}),$$

$$\hat{f}(y)={\hat{\mathbf{\Theta }}}^T\mathbf{\Phi }(\text{y})+{\mathbf{\epsilon }}_f,$$

$$\mathbf{\Phi }(\text{y})=\{\text{y}(k-d)$$

$$\text{u}(k-d)\cdots \text{u}(k-d-5)\},$$

$$\dot{\text{y}}$$

$$\text{A}+\text{B}\mathbf{\Lambda }{\text{K}}_y^T$$

$$\text{A}+\text{B}\mathbf{\Lambda }{\hat{\text{K}}}_y^T-{\text{A}}_m$$

$$\text{B}\mathbf{\Lambda }{\hat{\text{K}}}_r^T-{\text{B}}_m$$

$$\dot{\text{e}}(k)={\text{A}}_m\text{e}(k)+\text{B}\mathbf{\Lambda }[{\tilde{\text{K}}}_r^T\text{r}+{\tilde{\text{K}}}_y^T\text{y}-{\mathbf{\epsilon }}_f]$$

$$\text{V}(\text{e},{\tilde{\text{K}}}_y,{\tilde{\text{K}}}_r)$$

$$+\text{tr}({\tilde{\text{K}}}_r^T{\Gamma }_r^{-1}{\tilde{\text{K}}}_r^T|\mathbf{\Lambda }|)$$

$$\text{P}{\text{A}}_m+{\text{A}}_m^T\text{P}=-\text{Q},$$

$$\dot{\text{V}}(\text{e},{\tilde{\text{K}}}_y,{\tilde{\text{K}}}_r)$$

$$+2\text{tr}({\tilde{\text{K}}}_r^T{\mathbf{\Gamma }}_r^{-1}{\dot{\hat{\text{K}}}}_r|\mathbf{\Lambda }|)$$

$$={\text{e}}^T(\text{P}{\text{A}}_m+{\text{A}}_m^T\text{P})\text{e}+2{\text{e}}^T\text{P}\text{B}\mathbf{\Lambda }\left({\tilde{\text{K}}}_y^T\text{y}+{\tilde{\text{K}}}_r^T\text{r}-{\mathbf{\epsilon }}_f(\text{y})\right)$$

$$+2\text{tr}\left({\tilde{\text{K}}}_y^T{{\mathbf{\Gamma }}_{\mathbf{y}}}^{-1}{\dot{\hat{\text{K}}}}_y|\mathbf{\Lambda }|\right)+2\text{tr}\left({\tilde{\text{K}}}_r^T{\mathbf{\Gamma }}_r^{-1}{\dot{\hat{\text{K}}}}_r|\mathbf{\Lambda }|\right)$$

$$\begin{array}{rl}& =-{\text{e}}^T\text{Qe}-2{\text{e}}^T\text{PB}\mathbf{\Lambda }{\mathbf{\epsilon }}_f(\text{y})+2{\text{e}}^T\text{PB}\mathbf{\Lambda }{\tilde{\text{K}}}_y^T\text{y}\\ & +2\text{tr}\left({\tilde{\text{K}}}_y^T{{\mathbf{\Gamma }}_{\mathbf{y}}}^{-1}{\dot{\hat{\text{K}}}}_y\right)+2{\text{e}}^T\text{PB}\mathbf{\Lambda }{\tilde{\text{K}}}_r^T\text{r}+2\text{tr}\left(\Delta {\text{K}}_r^T{\mathbf{\Gamma }}_r^{-1}{\dot{\hat{\text{K}}}}_r\right)\end{array}$$

$$\begin{array}{rl}\dot{\text{V}}(\cdot )& =-{\text{e}}^T\text{Qe}-2{\text{e}}^T\text{PB}\mathbf{\Lambda }{\mathbf{\epsilon }}_f\\ & +2\text{tr}\left({\tilde{\text{K}}}_y^T({{\mathbf{\Gamma }}_{\mathbf{y}}}^{-1}{\dot{\hat{\text{K}}}}_y+\text{y}{\text{e}}^T\text{PB}sgn(\mathbf{\Lambda })\right)|\mathbf{\Lambda }|\\ & +2\text{tr}\left({\tilde{\text{K}}}_r^T({\mathbf{\Gamma }}_r^{-1}{\dot{\hat{\text{K}}}}_r+\text{r}{\text{e}}^T\text{PB}sgn(\mathbf{\Lambda })\right)|\mathbf{\Lambda }|\end{array}$$

$${\dot{\hat{\text{K}}}}_y=-{\mathbf{\Gamma }}_{\mathbf{y}}\text{y}{\text{e}}^T\text{P}\text{B}\text{sgn}(\mathbf{\Lambda }),{\dot{\hat{\text{K}}}}_r=-{\mathbf{\Gamma }}_r\text{r}{\text{e}}^T\text{P}\text{B}\text{sgn}(\mathbf{\Lambda }).$$

$$\dot{\text{V}}(\cdot )$$

$$\le -{\lambda }_{low}\Vert \text{e}{\Vert }^2+2\Vert \text{e}\Vert \Vert \text{PB}\Vert {\lambda }_{high}(\mathbf{\Lambda }){\mathbf{\epsilon }}_{max},$$

$$\chi =\left(\text{e}:\Vert \text{e}\Vert \le \frac{2\Vert \text{PB}\Vert {\lambda }_{high}(\mathbf{\Lambda }){\epsilon }_{max}}{{\lambda }_{low}(\text{Q})}\right).$$

$$f(\text{y}(k))=\text{arg }\underset{w}{\min }{\text{V}}_N(w,\mathbf{\Phi }(\text{y}))$$

$${\text{w}}_{k+1}\leftarrow \eta {\text{w}}_k-\alpha \sum _{i=1}^n{\nabla }_{\text{w}}\text{V}({\text{y}}_i,{\hat{\text{y}}}_i({\theta }_k)).$$

$${\text{y}}_m(t)=e^{{\text{A}}_{m}t}{\text{y}}_m(0)+{\int }_0^t{e}^{{\text{A}}_m(t-\tau )}{\text{B}}_m\text{r}(\tau )d\tau .$$

$$\begin{array}{r}{\text{A}}_m=\left[\begin{array}{ccc}-\frac{1334}{1705}& 0& 0\\ 0& -\frac{1334}{1705}& 0\\ 0& 0& -\frac{1334}{1705}\end{array}\right].\end{array}$$

Full text

I Introduction

Radiation-based treatment of head and neck (H&N) cancers often involve intensity-modulated radiotherapy (IMRT), which modulate radiation dosage and shaping of treatment beam to the precise size of tumor cells. IMRT carefully targets organs, while minimizing toxicity and exposure of organs at risk. Used alongside image guided radiotherapy (IGRT), IMRT assures precision of dosage targets: a patient is positioned on a treatment table after dosage planning by a physician, then a rigid robotic couch is used for motion alignment during surgery. While conventional RT uses rigid immobilization techniques such as masks, frames, arm positioning devices or vacuum mattresses [1], IGRT methods employ ultrasound, 3D imaging systems, 2D X-ray devices and/or computed tomography to instantly amend positioning errors, and improve daily radiotherapy fractions’ precision.

Precise, and repeatable patient positioning is therefore crucial in RT treatments when escalation of dose is necessary in a target volume and exposure of adjacent organs is to be minimized. Because it does not require rigid masks and body fixators, IGRT is more comfortable for the patient as well as more accurate with the aid of highly accurate localization systems. Sterzing [2], after an extensive review of IGRT methods, found IGRT to be more precise and safer than conventional radiotherapy. But setup errors (interfractional) or patient motion (intrafractional) errors often need to be accounted for during RT. While intrafractional errors can be minimized by highlighting the importance of voluntary stillness to the patient, suitable means of immobilization and adaptive positioning are necessary when the patient moves involuntarily or sleeps. This would assure precise and accurate targeting of critical organs while keeping the patient comfortable during treatment. Frameless and Maskless (F&M) IGRT is promising because it minimizes invasiveness and reduces setup times while comfortably positioning the patient.

We present Soft-NeuroAdapt, a set of three soft actuators that employ a neuro-controller to adaptively compensate for positioning deviations. Soft robot systems are elastomeric active and reactive nonlinear compliant systems, suitable to the task of human-machine interaction. They can dynamically change their stiffness and activate their surfaces to provide a desired motion on a human body part. Their characteristic morphological computation properties [3] make them suitable for autonomous control including inflation, deflation and reaching. Exhibiting highly nonlinear dynamics, controlling them for precise actuation is complicated. While a metric-based modeling approach may work in highly-structured environments, system parameters and dynamics change with different patients’ head and upper torso anatomy such that achieving precise control is best achieved with a learning-based method. But reliable high parameter-space adaptive-control methods rely on open-loop adaptation mechanisms based on extensive design techniques. Our goal is to derive a learning-based controller, going beyond task-specific, expert-driven methods in order to adaptively generalize to new control systems.

Following our recent investigative studies, [4, 5], on 1-DoF soft-robot compensation systems, we present a 3-DoF soft robot system that addresses 3-DoF involuntary intrafractional motions of the head and neck (H&N) region during F&M RT. This work presents a neuro-dynamics estimator that learns the underlying system model, then adapts to the model in its control law to provide bounded tracking of set trajectory. The controller corrects intrafractional (involuntary) patient motions along three defined axes namely head pitch, roll and elevation angles for a patient lying in a supine position on a table. We use three inflatable air bladders (IABs), actuated through a system of inlet and outlet solenoid proportional valves. The controller uses state feedback to provide bounded stability of states, a reference trajectory component to provide command tracking and a neural-network component to adaptively converge states that start outside of the sphere of stability into the region of stability. We perform head pose tracking in the 3D space of a stereo-camera, and we conduct experiments to validate the proposed bio-pneumatic system and controller.

This paper is organized as follows: section II presents related work, section III describes the hardware setup, section IV describes the vision segmentation algorithm and 3-DOF pose estimation, while section V describes the adaptive-neuro control algorithm. We discuss experiments in VI and provide discussions and future work in VII.

II Related Work

Studies on the feasibility of non-rigid immobilization devices during cranial stereotactic radiosurgery (SRS) fall largely under inspection-based non-rigid immobilization devices, and correction-based non-rigid immobilization systems.

Inspection-based Systems: In [6], Cervino et al. used expandable foams that fit a patient’s head while leaving the face opened. Patient set-up was performed using computed tomographic (CT) scans before treatment. They found an average treatment time of 26 minutes with patients who slept during experiments taking longer as a result of involuntary movements. In [7], the authors evaluated the accuracy of a head mold that minimally immobilized a patient’s H&N region while leaving the face free in a controlled positioning experiment with volunteers. 3D surface reconstruction imaging system was used in monitoring patients’ position so that treatment was stopped whenever motion exceeded a defined threshold. While the monitoring system showed great clinical accuracy, it required high patient cooperation to achieve immobilization. Using a head mold, an open face-mask, and a mouthpiece, Li et al. [8] quantified the residual rotation and positioning errors in an open-loop setting to ascertain the reduction in setup time during patient positioning setup. They reported the head mold and open face mask system could restrict head motions to within $0.6^{\circ}\pm 0.3^{\circ}$ with the time spent on in-situ motion corrections limited to $2.7\pm 1.0$ min.

Correction-based Systems: In [9], a 4-DOF motion correction system was used in a real-time motion compensation of a phantom and in human trials. An optical sensing system tracked the pose of the head and a decoupling control law regulated the $xyz$-translational and pitch motions of the head. They reported target accuracy of 0.5mm and $0.2^{\circ}$ with the decoupling controller for each axis of the 4-DOF motion. The system relied on stepper motor actuators and other electro-mechanical (EM) components positioned under the patient’s H&N region. These devices have the undesirable effect of attenuating treatment beams during RT, and as a result are not recommended during clinical procedure. The presence of the EM stages can significantly reduce the efficiency of the incident radiation targeted to tumor cells.

To address these concerns, we proposed soft-robot actuators for H&N motion compensation during F&M cancer radio-therapy [4, 5] actuated by pneumatic valves which are kept away from the head. Plastic hoses and silicone connectors convey air into and out of the soft robots that laid beneath the head so that our motion correction system does not interfere with planned dosage from the treatment beam. In our evaluations, we used pneumatic valves mounted on a movable plywood and connect these to the soft robots with long silicone tubes so that the medical physicist can separate EM components from the patient as deemed necessary during treatment. While our previous works related to controlling the 1-DoF motion of the head along the z-axis, here we address the 3-DOF control of the head using a surface reconstructing sensor and control the pose of the head to follow set trajectories in real-time.

III Hardware Design Overview

The actuation mechanism consists of three custom-designed inflatable air bladders (IABs) made from elastomeric polymers. The base IAB (beneath the head’s posterior) is 180mmx280mm when flat and inflates to a maximum width of $\sim 75$ mm while the other two are 180mmx140mm in size. The IABs consist of inflatable rubber encased in a breathable foam pad for comfort, modified to be the size of an average adult male; they have separate inlet and outlets openings, connected with crack-resistant polyethylene tubing (1/8” ID and 1/4” OD); this sustains pressure of up to 320psi.

Each hose leads to a proportional solenoid valve, which is in turn connected to rectangular manifolds (one manifold to the inlet supply and the second to the outlet supply). We use six Dakota Instruments EM valves (Model PSV0105, Orangeburg, NY, USA) to supply proportional torques to the soft actuators. A regulated air canister supplied constant air pressure at 15psi to the inlet-air conveying manifold, while a suction pump supplied vacuum pressure at 12 psi to the valves that removed air from the bladders. The air rate of flow into or out of each bladder was controlled via custom-built voltage regulating circuits which got PWM signals from a National Instruments (NI) myRIO microcontroller.

We 3D printed a custom manikin head, measuring $155\times 240\times 200$ mm ($W\times L\times D$) and comparing between 50% and 75% weight of a typical adult male head or 99% of a typical adult female head; it weights 5kg – above average and reasonable. The head was fitted with a ball-joint in the neck to replicate motion of the human head about the neck. An Ensenso 3D camera is mounted near $45^{\circ}$ above the head to measure the pose of the head in real time. All vision processing, systems modeling and control laws were computed on a CORSAIR PC. We exchange the neuro-control and sensor signals via the publish-subscribe IPC of the ROS middleware installed on the PC. Adaptive control laws were sent via udp packets to the RIO microcontroller. The system setup is shown in Fig. 1.

The reference frame of the head is described as follows: the pitch/x-axes points from the left ear out of the right ear, the yaw/z-axes points from the back of the head through the forehead through, and the roll/y-axes goes from the neck through the top of the head. The left and right bladders control the roll angles/x-axes motions while the bladder underneath the head, henceforth referred to as the base bladder, controls the pitch angles and z-axis motions. The reference frame is illustrated in Fig. 2.

IV Vision-based Pose Estimation

The model head lies in a supine position above a planar table as shown in Fig. 1. We employed a 3D camera from Ensenso GmbH (model N35) to reconstruct the surface image and measure head pose. The N35 camera captures multiple image pairs during exposure; each image pair is made up of different patterns, controlled by piezo-actuators. A stereo-matching algorithm gathers the information from all image pairs after capture to produce a high-resolution point cloud (PCL) of the scene [10]. We mounted the 3D sensor such that its lens faced the head at approximately $45^{\circ}$ from the vertical during experiments. Our goal is to control the motion of the head about three axes, namely $z$, pitch and roll axes as described in section III; this section presents how we go about extracting representable features from the face.

IV-A Face Segmentation

The dense point cloud of the scene has (i) marked jump in rendered points along the z-axis of the camera; this is because of the single view angle by the camera; (ii) the scene clutter and lack of multiple camera view angles does not affect the representation of the face; (iii) thus, through spatial decomposition of the scene, we can separate the face from the scene. However, the point cloud is computed from monochromatic IR image pairs (with no texture information) making morphological operations difficult; due to the multiple image pairs used in 3D reconstruction to generate a highly accurate measurement, the camera is limited to a maximum frame rate of 10Hz. Inspired by Rusu’s work [11], we divide the segmentation problem into stages, with each stage involving segmenting out candidates that do not belong to the object we want to identify (the frontal face) in the scene. Our engineering philosophy in the segmentation phase is inspired by spatial decomposition methods that determine subdivisions and boundaries to allow retrieval of data that we want given a measure of proximity. In this case, we know that the location of the table cannot exceed a given height during experiments and the camera’s position is fixed while the head moves based on bladders’ actuation. Separating objects that represent planar 2D geometric shapes from the scene therefore simplifies the face segmentation algorithm. By finding and removing objects that fit primitive geometric shapes from the scene, clustering of the remaining objects would yield the face of the patient in the scene. We fit a simplified 2D planar object to the scene such that searching for points $\textbf{p}_{i}\in\mathcal{P}$ that support a 2D plane can be found within a tolerance defined by the inequality $0\leq|d|\leq|d_{max}$, where $|d_{max}|$ represents a user-defined threshold to segment out [11].

We proceed as follows:

• •

The point cloud of the scene was acquired from the computed disparity map of the two raw camera images;

• •

To minimize sensor noise whilst preserving 3D representation, the acquired point cloud was downsampled using a SAmple Consensus (SAC)-based robust moving least squares algorithm (RMLS) [11, §6];

• •

We then searched for the edges of 2D planar regions in the scene with Maximum Likelihood SAmple Consensus (MLESAC)[12], and we bound the resulting plane indices by computing their 2D convex hull;

• •

A model fitting stage extrudes the computed hull (of objects lying above the 2D planar region) into a prism model based on a defined $L_{1}$ Manhattan distance; this gives the points whose height threshold is about the region of the face in the scene [13];

• •

We then cluster the remaining points based on a heuristically determined $L_{2}$ distance between points remaining within the polygonal plane. The largest cluster gives us the face.

The result of the resampling algorithm is shown in the top-right image of Fig. 3.

To simplify the complexity of the planar structure in the scene, the table is modeled as a 2D planar geometric primitive so that finding points that fit a defined model hypothesis involves estimating a single distance to the frontal plane of the table surface rather than multiple points if the model was represented with points. Searching for horizontal planes that are perpendicular to the z-axis of the head is carried out using the maximum likelihood SAC [12] algorithm implemented in the PCL Library [14] to generate model hypotheses. The plane segmentation algorithm is defined in Table I.

The plane segmentation process is run once. Once the plane model is found, its indices and those of objects lying above it are separated and stored in separate data structures. Every subsequent iteration consists of (i) computing the 2D convex hull of point indices of objects above the table using the Qhull library111The Qhull library: http://www.qhull.org/, (ii) using a pre-defined prismatic model candidate to hold extruded points to the approximate facial height above the table; and (iii) separating the face from every other point in the resulting cloud through the Euclidean clustering (EC) method of [15]. A distinct point cluster is defined if the points in cluster $C_{i}=\{\textbf{p}_{i}\in\mathcal{P}\}$ and cluster $C_{j}=\{\textbf{p}_{j}\in\mathcal{P}\}$ satisfy the $L_{2}$-distance threshold

[TABLE]

Finding the face in the scene after carrying out EC algorithm is a question of finding the largest index in the list $\mathcal{C}$. This takes $O(n)$ (linear) time for $n$ clusters. The face segmentation results are presented in Fig. 3. We then compute the Cartesian position of the face with respect to the camera origin by taking the center of mass of the segmented facial region (bottom-right image of Fig. 3). This is obtained by calculating the mean-value of all the points in the resulting cloud ($\approx$ 600 points on average).

IV-B Head Pose Estimation

With the facial point cloud segmented, we define three points on the head. Our goal is to compute the optimal translation and rotation of the head from a model point set $\textbf{X}=\{\overrightarrow{x}_{i}\}$ to a measured point set $\textbf{P}=\{\overrightarrow{p}_{i}\}$, where $N_{x}=N_{p}=3$, and the point $\overrightarrow{x}_{i}\in\textbf{X}$ has the same index as $\overrightarrow{p}_{i}\in\textbf{P}$. Following the approach of Besl and McKay in [16], we compute the cross-covariance matrix of P and X as $\Sigma_{px}$, extract the cyclic components of this skew symmetric matrix as $\bm{\Delta}$, and use it to form the symmetric $4\times 4$ matrix $\textbf{Q}(\bm{\Sigma}_{px})$ as follows,

[TABLE]

The unit eigenvector, $q_{R}$, that corresponds to the maximum eigenvalue of $\textbf{Q}(\bm{\Sigma}_{px})$ is selected as the optimal rotation quaternion; we find the optimal translation vector as

[TABLE]

where $\mu_{x}$ and $\mu_{p}$ are the mean of point sets X and P respectively. Obtaining the roll, pitch and yaw angles from $q_{R}$ is trivial and the pose of the face is described by tuples $[q_{T},q_{R}]=\{x,y,z,\theta,\phi,\psi\}$ with respect to the world frame. Given the 3-DOF setup, we choose to control three states of the head: $z$, $\theta,\,\phi$ ($i.e.$ z, roll, and pitch).

V Control Design

Our control philosophy is governed by the state feedback and feedforward regulation problem, with an adaptation mechanism based off an estimation of the head pose given a priori information about the system’s states and past control actions. We propose an adaptive control strategy in a Bayesian setting, which given an initial prior distribution of controls and 3-DOF head pose, minimizes a cost criterion as the expected value of control laws that will yield a future desired head pose. We consider the pwm voltages that power the valves as input, u, the head pose as the output, y and an unknown disturbance $\textbf{w}(k)$. We first describe the nonlinear function approximator model $\hat{f}\left(\textbf{u}(k-d),\textbf{y}(k),\textbf{w}(k)\right)$, which is constructed from memory-based input and output experimental data that satisfy

[TABLE]

that satisfy the Lipschitz continuity. (3) implies an input $\textbf{u}(\cdot)$ at time $k-d$, produces an output $\textbf{y}(k)$ at $d$ time instants later. The next section describes how we formulate the class of minimum error variance controllers that predict the effect of actions $\textbf{u}(\cdot)$ on states $\textbf{y}(\cdot)$ using a self-tuning regulator.

V-A Adaptive Neuro-Control Formulation

Following our previous approach in [5, §IV.B], we fix a persistently exciting input signal $u_{ex}\in L_{2}\cap L_{\infty}$ to excite the nonlinear modes of the system. We then parameterized the system with a neural network with sufficient number of neurons. The neural network (NN) provided information on the changing parameters of the system during control trials. The adjustment mechanism is computed from inverse Lyapunov analysis, where we choose adaptive laws that guarantee a nonpositive-definite Lyapunov function candidate when evaluated along the trajectories of the error dynamics.

Our contribution is the approximation of the nonlinear system by a long short-term memory (LSTM) [17], equipped with an adequate number of neurons in its hidden layers. We parameterized the last layer of the network with a fully connected layer that outputs control torques to the valves. The neural network can be seen as a memory-based model that remembers effective controls for the adaptation mechanism in the presence of uncertainties and external disturbance.

The neural network is shown in Fig. 4. Depending on the region of attraction of the system the network is approximating, it parameterizes the nonlinear dynamical system $f(\cdot)$ and maps the parameterized model to appropriate valve torques. There exists additional feedforward + feedback terms in the global controller (introduced shortly) that guarantee system stability and robustness to uncertainties. Therefore, the global controller keeps the states of the system bounded under closed-loop dynamics, ensures convergence to desired trajectories from states that are initialized outside the domain of attraction, and guarantees robust reference tracking in the presence of non-parametric uncertainties.

For the multi-input, multi-output (MIMO) adjustable system,

[TABLE]

where $\textbf{y}\in R^{n},\textbf{u}\in R^{m}$ are known input and output vectors, and $\textbf{A}\in R^{n\times n},\bm{\Lambda}\in R^{m\times m}$ are unknown matrices, $\textbf{B}\in R^{n\times m}\,,sgn(\bm{\Lambda})$ are known matrices, and $\textbf{w}(k)\in R^{n}$ is a bounded time-varying unknown disturbance, upper-bounded by a fixed positive scalar $\textbf{w}_{max}$. We make the following assumptions:

• •

a dynamic RNN with $N$ neurons, $\varphi(\textbf{y})$, exists that will map from a compact input space $\mathbb{U}\supset\textbf{u}$ to an output space $\textbf{y}\subset\mathbb{Y}$ on the Lebesgue integrable functions with closed interval $[0,T]$ or open-ended interval $[0,\infty)$;

• •

the nonlinear function $f(\cdot)$ is exactly $\Theta^{T}\Phi(\textbf{y})$ with vectorized coefficients $\bm{\Theta}\in R^{N\times m}$ and a Lipschitz-continuous vector of basis functions $\bm{\Phi}(\textbf{y})\in R^{N}$;

• •

inside a ball $\textbf{B}_{R}$ of known, finite radius $R$, the ideal NN approximation $f(\textbf{y}):R^{n}\rightarrow R^{m}$, is realized to a sufficient degree of accuracy, $\bm{\varepsilon}_{f}>0$;

• •

the process noise $\textbf{w}(k)$ is estimated alongside model parameters by the dynamic RNN;

• •

outside $\textbf{B}_{R}$, the NN approximation error can be upper-bounded by a known scalar function $\bm{\varepsilon}_{max}$ such that $\|\varepsilon\|\leq\varepsilon_{max},\quad\,\forall\quad\textbf{y}\in\textbf{B}_{R}$;

• •

there exists an exponentially stable reference model

[TABLE]

with a Hurwitz matrix $\textbf{A}_{m}\in R^{n\times n}\text{ and }\textbf{B}_{m}\in R^{n\times m}$ commanded by a reference signal $\textbf{r}\in R^{m}$. For this system, we note that $n=3$ and $m=6$. Our objective is to design an model-reference adaptive controller (MRAC) capable of operating in the presence of parametric ($\bm{\varepsilon}_{f}$), and non-parametric ($\textbf{w}(k)$) uncertainties so as to assure the boundedness of all signals within the closed-loop system. We propose the following controller

[TABLE]

where $\hat{\textbf{K}}_{y}$ and $\hat{\textbf{K}}_{r}$ are adaptive gains to be designed shortly. The $\hat{\textbf{K}}_{y}^{\textbf{T}}\textbf{y}$ term keeps the states of the approximation set $\textbf{y}\in\textbf{B}_{R}$ stable, while the $\textbf{K}_{r}^{\textbf{T}}\textbf{r}$ term causes the states to follow a given reference trajectory. The function approximator $\hat{f}(\cdot)$ ensures states that start outside the approximation set $\textbf{y}\in\textbf{B}_{R}$ converge to $\textbf{B}_{R}$ in finite time (it converges non-parametric errors $\bm{\varepsilon}_{f}$ that puts certain states out of the approximation set into $\textbf{B}_{R}$). We can generally write the NN model as

[TABLE]

where $\hat{\bm{\Theta}}^{T}$ denotes the vectorized weights of the neural network and $\bm{\Phi}(\textbf{y})$ denotes the vector of inputs and outputs defined as

[TABLE]

and $\bm{\varepsilon}_{f}$ is the approximation error. The closed-loop dynamics therefore become

[TABLE]

We assume nonlinear function and approximator matching conditions, $f(\cdot)=\hat{f}(\cdot)$, such that after rearrangement, (8) can be written as

[TABLE]

Furthermore, we assume model matching conditions with ideal constant gains $\textbf{K}_{y}\text{ and }\textbf{K}_{r}$ so that

[TABLE]

from which

[TABLE]

where $\tilde{\textbf{K}}_{y}^{T}=\textbf{K}_{y}^{{}^{T}}-\hat{\textbf{K}}_{y}^{T}$ and $\tilde{\textbf{K}}_{r}^{T}=\textbf{K}_{r}^{{}^{T}}-\hat{\textbf{K}}_{r}^{T}$. The generalized error state vector $\textbf{e}(k)=\textbf{y}(k)-\textbf{y}_{m}(k)$ has dynamics $\dot{\textbf{e}}(k)=\dot{\textbf{y}}(k)-\dot{\textbf{y}}_{m}(k)$, so that by substituting (5) and (8) into $\dot{\textbf{e}}$, we have

[TABLE]

The estimation error will be bounded as long as $\textbf{y}\in\textbf{B}_{R}$. Our goal is to keep $\textbf{y}\in\textbf{B}_{R}$.

Theorem: Given correct choice of adaptive gains $\hat{\textbf{K}}_{y}$ and $\hat{\textbf{K}}_{r}$, the error vector $\textbf{e}(k)$, with closed loop time derivative given by (12) will be uniformly ultimately bounded, and the state y will converge to a neighborhood of r.

Proof: We choose a Lyapunov function candidate V in terms of the generalized error state space e, gains, $\tilde{\textbf{K}}_{y}^{T},\,\tilde{\textbf{K}}_{r}^{T}$, and parameter error $\bm{\varepsilon}_{f}(\textbf{y}(k))$ space ([18], [19], [20]) as follows

[TABLE]

where $\bm{\Gamma_{y}}$ and $\bm{\Gamma}_{r}$ are fixed symmetric, positive definite (SPD) matrices of adaptation rates, $\textbf{tr}(\textbf{A})$ denote the trace of matrix A and P is a unique SPD matrix solution of the algebraic Lyapunov function

[TABLE]

where Q is a SPD matrix. Take the time derivative of (13)

[TABLE]

[TABLE]

[TABLE]

Since $x^{T}\,y=\textbf{tr}\left(y\,x^{T}\right)$ from trace identity, we have

[TABLE]

where for a real-valued $x$, we have ${x}=\text{sgn}(x)|x|$. The first two terms in (15) will be negative definite for all $\textbf{e}\neq 0$ since $\textbf{A}_{m}$ is Hurwitz and the other terms in (15) will be identically null if we choose the adaptation laws

[TABLE]

The time-derivative of the Lyapunov function can then be written as

[TABLE]

where $\lambda_{low},\lambda_{high}$ represent the minimum and maximum characteristic roots of $Q\text{ and }\bm{\Lambda}$ respectively. $\dot{\textbf{V}}(\cdot)$ is thus negative definite outside the compact set

[TABLE]

and we conclude that the error e is uniformly ultimately bounded. As e converges to a neighborhood of 0, y converges to a neighborhood of $\textbf{y}_{m}$. From the stable model reference system in (5), y converges to a neighborhood of r. Note that asymptotic convergence of e to zero is not guaranteed but the parametric errors are guaranteed to stay bounded.

V-B Network Design

We require accurate mapping of temporally lagged patterns in inputs to output states, a dynamic nonlinear model of valve encoder values to sensor measurements that accurately maps $f(\cdot)$ in (4). We choose a LSTM [17] due to its capacity for long-term context memorization and inherent multiplicative units that avoid oscillating weights or vanishing gradients when error signals are backpropagated in time [21, 17]. LSTMs truncate gradients in the network where it is harmless by enforcing constant error flows through their constant error carousels. As a result, LSTMs are robustly more powerful for adaptive sequence-to-sequence modeling or mapping data that temporally evolve in time. Their biological model makes them more suitable for adaptive robotics such as soft robots than previously used artificial NNs such as feedforward networks [22], radial basis-functions [23, 20] or vanilla RNNs [24].

The NN model takes a memory-based concatenated vector of current inputs and past outputs as in (7), propagates them through three hidden layers, with each layer made up of $\{9,6,6\}$ neurons each, applies 30% dropout and then maps the last layer to a fully connected layer that generates valve torques. The architecture of the neuro-controller is shown in Fig. 4. The last layer is designed to generate appropriate valve torques based on an internal model of the plant. A self-tuning adaptive control law (with a feedforward regulation and state feedback component) adapts to the internal parameters of the plant to ensure stability of the system and bounded tracking of given trajectory. The overall network has neuron connection weights and thresholds of approximately 1,400. This makes search for a suitable controller feasible.

The LSTM model estimates a model $f(\textbf{y})$, that minimizes the mean-squared error between predicted output $\hat{\textbf{y}}(k)$ and actual output $\textbf{y}(k)$ according to

[TABLE]

$\text{ where }\textbf{V}_{N}(\textbf{w},\bm{\Phi}(\textbf{y}))=\sum_{t=1}^{\mathcal{K}}\sum_{i=1}^{n}\dfrac{1}{2}(\hat{\textbf{y}}_{i}(t)-\textbf{y}_{i}(t))^{2},$ and $\Phi(\textbf{y})$ is a regression vector as defined in (7) on a bounded interval $[1,N]$. (19) is minimized using stochastic gradient descent so that at each iteration, we update the parameters (weights) of the network $\textbf{w}_{i}$ based on the ordered derivatives of $\textbf{V}_{N}(\textbf{w},\bm{\Phi}(\textbf{y}))$ (Werbos [25]) $i.e.$

[TABLE]

$\eta$ (set to 1) hastens the optimization in a direction of low but steepest descent in training error, and $\alpha$ is a sufficiently small learning rate (set to $5\times 10^{-3}$), and $\nabla_{w}\textbf{V}(\theta,\bm{\Phi}(\textbf{y}))$ is the derivative of V with respect to $w$ averaged over the $k$-th batch (we used a batch size of 50). We initialized the weights of Fig. 4 from a one-dimensional normal distribution with zero-mean and unit variance.

VI Experiment: Soft actuators Control

The headpose is determined based on our formulation in IV. The 3-DOF pose of the head is made up of the state tuple $\{z(k),\theta(k),\phi(k)\}$.

VI-A Adaptive Control Parameters

We sample from the parameters of the trained network and we set $\hat{f}(\cdot)$ in (6) to the fully connected layer of samples from the network. We publish the control law from the neural network and subscribe in a separate node. The gains $\hat{\textbf{K}}_{y}$ and $\hat{\textbf{K}}_{r}$ in (16), were found by solving the ODEs iteratively using a single step of the integral of the solutions to $\dot{\hat{\textbf{K}}}_{y}(t),\,\dot{\hat{\textbf{K}}}_{r}(t)$. Our solution is an implementation of the Runge-Kutta Dormand-Prince 5 ODE-solver available in the Boost C++ Libraries222https://goo.gl/l7JyYe. We found a step-size of $0.01$ to be realistic. $\textbf{y}_{m}$ in (5) is computed based on the solution to the forced response of the linear system,

[TABLE]

We set $\textbf{y}_{m}(0)=\textbf{y}(0)$ at t = 0 and for a settling time requirement of $T_{s}=5secs$ at which the response remains within $2\%$ of final value, we find that

[TABLE]

For a nonnegative Q and a positive definite P, the pair $(\textbf{Q},\textbf{A}_{m})$ will be observable (LaSalle’s theorm) so that the dynamical system is globally asymptotically stable. After searching, we picked a positive definite $\textbf{Q}=diag(100,100,100)$ for the dissipation energy in (17) and set $\bm{\Lambda}=I_{3\times 3}$ so that solving the general form of the lyapunov equation, we have

[TABLE]

The six solenoid valves operate in pairs so that two valves create a difference in air mass within each IAB at any given time. Therefore, we set

[TABLE]

The non-zero terms in (23) denote the maximum duty-cycle that can be applied to the Dakota valves based on the software configuration of the NI RIO PWM generator.

VI-B Results and Analysis

The three DoFs of the head are coupled and there is a limited reachable space with the IABs. It is therefore paramount that desired trajectories be ascertained as physically realizable before rolling out control trials. We therefore placed the head to physically realizable positions in open-loop control settings before testing the close-loop control system on such feasible goal poses.

Fig. 5 show the performance of the controller when commanded to move the head from $[z,\theta,\phi]^{T}=[2.5mm,.25^{o},35^{o}]^{T}$ to $[14mm,1.6^{o},45^{o}]^{T}$. We observe strong steady-state convergence along 2-DoFs, namely z and pitch axes with a 20 second rise time. The roll motion is however characterized by offshoots that may be caused by the coupled DOF. We perform a second experiment, seen in Fig. 6, where we evaluate the performance of the controller on the roll angle of the head. We observe that the controller behaves well controlling the roll motion in isolation. The overshoot of Fig. 5 are likely due to coupled dynamics not accounted for in our formulation.

While our results are promising, showing the feasibility of the control law along the physically realizable axes of motions, the coupled degrees of freedom need further investigation to achieve independent and precise axial motion control whilst preserving global head goal requirements.

VII Discussion and Future Work

We have presented a soft robot motion compensation system that uses a robust adaptive neurocontroller to correct patient positioning deviation in F&M RT along 3 degrees of freedom. Unlike related works that perform an inspection-based correction mechanism or employ radiation attenuating devices such as stepper motors positioned underneath the patient’s head on the table, our method offers the elimination of reduction in the efficacy of dosage plans through compliant nonlinear soft elastomeric polymers. It also eliminates the pain of wearing metallic rings that rigidly immobilize a patient during treatment setups and dosage administration. The prospect of autonomously adapting to changing model parameters in our controller by learning compact and portable state representations of complex environments has widespread implications for autonomous robots.

Further research in this direction will focus on decoupling the control laws for the coupled states of the system whilst preserving global pose objectives. We are investigating control of custom-made multi-chambered air-bladders. Extensible to high-space DoF control, we will evaluate multiple soft-robot systems on a full-fledged human phantom and volunteer human trials.