Nonlinear Model Predictive Control of a Robotic Soft Esophagus

Strictures caused by esophageal cancer can narrow down the esophageal lumen, leading to dysphagia. Palliation of dysphagia has driven the development of a robotic soft esophagus (RoSE), which provides a novel <italic>in vitro</italic> platform for esophageal stent testing and food viscosity studies. In RoSE, peristaltic wave generation and control were done in an open-loop manner since the conduit lacked visibility and embedded sensing capability. Hence, in this work, RoSE version 2.0 (RoSEv2.0) is designed with embedded time of flight (TOF) and pressure sensors to measure conduit displacement and air pressure, respectively, for modeling and control. Model predictive control (MPC) of RoSEv2.0 is implemented to govern the peristalsis and air pressure profile autonomously. The implemented MPC used sparse identification of nonlinear dynamics with control (SINDYC) models to estimate the future states of ROSEv2.0. The dynamic models are discovered from the TOF and pressure sensor data. Peristalsis waves of speed 20 mm <inline-formula><tex-math notation="LaTeX">$\cdot$</tex-math></inline-formula> s<inline-formula><tex-math notation="LaTeX">$ ^{-1}$</tex-math></inline-formula>, wavelength 75 mm, and amplitudes 5, 7.5, and 10 mm were successfully generated by the MPC. Additionally, RoSEv2.0 with the MPC was employed to perform stent migration testing with various food boluses consistencies. The major contribution claimed in this article is the application of SINDYC-based MPC to solve the closed-loop control problem of RoSE for achieving desired peristaltic waves.

peristaltic waves is a significant shortcoming of concern to the patients [3].
The migration can be minimized by improving stent designs, but evidence to determine which stent design is better than the other is confined to stent analytical [4] and numerical modeling [5] studies because in vivo stent testing can raise major ethical concerns. Garbey et al. [5] and Mozafari et al. [6] are two relevant numerical studies that have performed stent testing on numerical esophageal models. While the former research work focused on stent migration, the latter identified stent characteristics that significantly impacted stent migration.
A robotic soft esophagus (RoSE) has been developed as an alternative platform for conducting stent testing under bolus swallow conditions [7]. Since RoSE has been designed to mimic the human swallow behavior, matching the attributes of the biological esophagus; thus, instead of actual patients, RoSE has been used to conduct the study on the various stent designs before implanting them in patients with malignant esophageal strictures [7].
The dynamic modeling of soft and continuum robots, such as the RoSE, has been challenging and remains an active research topic [8]. For soft continuum robots' controller design, many sensors are typically required to measure the high-dimensional states to provide the controller with state feedback. Hence, several authors have considered task-space closed-loop control that does not require a dynamic model and a large number of sensors and can provide a robust and simplistic solution to the control problem without any previous information of the robots' dynamics. While Lee et al. [9] proposed a model-less control framework for a soft fluid-driven continuum robot to follow a defined path accurately, Li et al. [10] introduced an adaptive Kalman-based controller to estimate the Jacobian of a pneumatic muscle actuator driven continuum robot.
To avail from the vast knowledge of prevailing automatic control theory, soft robotics dynamic modeling is required since the more information a controller has about the plant, the better the tracking result [11]. Additionally, unlike model-free control, model-based control approaches also take the compliant nature of such robots into consideration, and the controller can be operated in a high-speed environment. There is a growing body of the literature that recognizes the importance of model-based control techniques, such as optimal control [12] and model predictive control (MPC) [13], because they are capable of routinely taking actual plant constraints into account in time.
MPC is an advanced and well-known technique of controlling highly nonlinear processes with constraints [13], [14]. Based on 0278-0046 © 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
the past observations, MPC employs the concept of receding horizon (forward shifting of prediction horizon) to anticipate potential disturbances and thereby to generate control sequence that is a proposition among possible alternatives that can manage disturbances well at the next time step [13], [15]. Unlike traditional control theory, MPC is effective if there is a process deadtime or if the setpoint/reference trajectory is well defined ahead of time. Additionally, it provides the flexibility to formulate and modify the objective function for optimal control. Since the 1980s, there is a dramatic increase in the application of MPC in industries [15], such as power systems [16] and robotics [17], [18], because MPC can model and thereby predict nonlinear and nonminimum phase dynamics systems very accurately. Advanced control techniques, such as the MPC, relies on developing a suitable plant model for optimum performance. Soft robots, such as the RoSE [19] and soft robotic gastric simulator [20], manufactured by soft elastomers, such as silicone rubber, cannot be modeled from the perspective of first principles due to their complex, continuous, and highly compliant intrinsic deformation, which leads to infinite degrees of freedom [21]. With the increasing complexity of the soft robots and rigid performance specifications, the researchers have revealed that MPC implementation with improved predictive model can significantly enhance the control performance of complex nonlinear systems [18], [22].
Data-driven techniques are not a panacea for all engineering problems, but if used correctly, it provides a way to leverage the collected data from experiments maximally. The sparse identification of nonlinear dynamics with control (SINDYC) approach [23], [24], has succeeded in addressing challenges in the data-driven modeling of RoSE [25]. Unlike traditional black-box data-driven modeling techniques, SINDYC has demonstrated the power to identify nonlinear dynamics from data captured from a nonlinear system. The method has been applied and validated with several dynamical systems, such as the damped harmonic oscillator with linear or cubic dynamics, the Lorenz system, and identifying bifurcations and normal forms [23]. Fields, such as soft robotics, where there is the absence of governing equations, may benefit from methods, such as the SINDYC. Bhattacharya et al. [25] has proven that the SINDYC algorithm promotes sparsity by identifying a few linear and nonlinear terms governing the dynamics of the RoSE. Additionally, unlike other machine learning (ML) techniques, such as the artificial neural networks, SINDYC demands less training data, less execution time, and avoids overfitting [26]. Integrating the RoSE SINDYC model in the MPC framework can take advantage of the underlying dynamics in achieving the prescribed peristalsis profile for stent testing.
The primary aim of this article is to actuate the RoSE conduit with predefined peristaltic wave trajectories by implementing SINDYC-based MPC framework. The fulfillment of the aim will make RoSE capable of generating various peristaltic wave shapes, which will further enhance the in vitro stent migration testing potential of RoSE. Five significant steps were followed to implement the MPC for stent migration testing. First, RoSE version 2.0 (RoSEv2.0) was developed with an embedded array of time-of-flight (TOF) sensors to address the issue of limited visibility inside the RoSEv2.0 conduit. Second, RoSEv2.0 was actuated with various pressure varying trajectories and data from TOFs and valve pressure sensors (VPSs) were collected. Third, discrete-time formulation of SINDYC referred to as discretetime SINDYC (DTSINDYC) was applied to the collected data to discover discrete time differential equation (DTDE) models of RoSEv2.0. Fourth, nonlinear MPC using the DTDE models were designed and implemented to control the peristaltic motion of the RoSEv2.0 conduit. Last, stent migration data for a candidate stent were recorded for various peristalsis speed and bolus swallow conditions.
The research work makes several noteworthy contributions to the field of soft robotics. This research work has solved the sensing issue for modeling and control of RoSE by introducing RoSEv2.0 (improved version of RoSE) with TOF sensors. Additionally, to the best of our knowledge, the SINDYC-based MPC framework presented in this research work has been applied for the first time in a soft robot. The closed-loop control in the form of MPC has been implemented for the first time in RoSEv2.0. Overall, the contribution claimed in this article is the application of SINDYC-MPC to solve the closed-loop control problem of RoSE for peristaltic wave tracking.
Additionally, in this research work, the presented modeling and control methodology is not a robot or sensor specific, and it can be extended to any soft-robotic system with any sensor provided input-output datasets can be generated. Hence, to verify the versatility of this approach, MPC of RoSEv2.0 controlling the chamber air pressure was implemented with VPS by following the same implementation methodology described with the TOF sensors.

II. PROBLEM FORMULATION
The foundation of supervised ML techniques, such as the SINDYC, is qualitative and quantitative data collection. To a greater extent, the quality of the training data determines the SINDYC modeling generalization when applied to any soft robot, such as RoSE [19]. The updated RoSE (RoSEv2.0) has 12 identical layers (L 1 to L 12 ) such that each layer contains four air pressure chambers surrounding the conduit axis symmetrically. Each horizontal whorl of four chambers at each layer is connected to an electropneumatic pressure valve (ITV-0030-3BS, SMC, Indianapolis, IN, USA) such that RoSEv2.0 has 12 independent inputs to control the actuation of RoSEv2.0. In contrast with RoSE, RoSEv2.0 is developed with embedded sensing capability. The significant variation between RoSE and RoSEv2.0 is the outer casing construction, which is replaced with a transparent liquid polydimethylsiloxane (PDMS) (SYL-GARD 182, Dow, Midland, MI, USA) casing in RoSEv2.0 (see Fig. 1). An array of TOF (VL6180X, STMicroelectronics, Geneva, Switzerland) sensors are placed on top of the PDMS layer to measure the conduit deformation laterally from the outside.
The DTSINDYC models M 1 and M 2 based on TOF and VPS, respectively (see Table I), for implementing the MPC are built on a small number of datasets collected from the RoSEv2.0 sensors  and aims to discover the underlying DTDEs given by (1) such that the equations best-fit the captured dataset In (1), u k ∈ R n and x k ∈ R n represent the digital values for pressure command applied to the pneumatic valves at kth time step to govern the RoSEv2.0 actuation and its corresponding TOF or VPS measurements, respectively. g(.) represents a nonlinear polynomial function consisting of first-and secondorder polynomial, and constant terms, mapping R n × R n → R n . Equation (1) also represents the predictive nature of the DTSINDYC models of RoSEv2.0, which will be used for implementing the MPC. At any time-instant k, given x k and u k of RoSEv2.0, the future of the RoSEv2.0 states denoted by x k+1 can be predicted. In this research work, g(.) represents two different DTSINDYC models M 1 and M 2 , generated from the VPS and TOF captured data, respectively. Table I provides the description and physical representation of the state variables associated with the models. The models are used in the MPC for controlling conduit displacement or chamber air pressure in case of peristaltic actuation [7]. Since at least three RoSEv2.0 layers are required to match the human esophagus peristalsis wavelengths [27], thus by emphasizing the TOF sensors measurement accuracy, layers L 5 , L 6 , and L 7 were chosen. The MPC of RoSEv2.0 using M 1 or M 2 is implemented to compute the control law u(.|x k ) = {u j+1 , . . . , u j+k , · · · , u j+N u } at any time-step j, for achieving a prescribed peristaltic actuation profile of the robot's conduit over a control horizon N u and prediction horizon N p , provided the current TOF or VPS measurement x j is available. Hence, the implicit feedback law for controlling the valves and, hence, controlling the RoSEv2.0 can be written as u(j + 1|x j ) = u j+1 (2) where u j+1 is the first row in the optimized actuation sequence, representing pressure commands applied to the L 5 , L 6 , and L 7 pneumatic valves, beginning at the initial condition x j . The entire RoSEv2.0 conduit L 1 to L 12 does not contain any rigid skeletal boundary, and when pressurized with air, the whorl of four chambers per layer inflate and occlude, mimicking circular muscle activation. The complete rigid-boundary less Ecoflex structure of the conduit is shown in Fig. 2(j). The pneumatic valves corresponding to each RoSEv2.0 layer are interfaced with a Raspberry Pi 4 Model B through an analog-to-digital converter and digital-to-analog converter (DAC) board. The pneumatic valves have a set pressure range of 1 − 500 kPa, operates with an input signal range of 0 − 10 V, characterizing a resolution of 50 kPa · V −1 . Since the 8-b DACs connected to the valves have a resolution of 0.02 V · step −1 thus, each digital value applied to the DAC-valve assembly will generate a pressure with 1 kPa · step −1 . For complete occlusion of the RoSEv2.0 conduit, maximum pressure of 47 kPa is required, which is regarded as the upper limit for the control signals generated by the MPC.
Custom firmware protocols for generating peristalsis, written in Python 3.7, are developed on the Pi to actuate the robot in open-loop, with continuous, independent variable pressure trajectories, as discussed in Section II. By utilizing the peristaltic protocols, training and validation data for the DTSINDYC models are collected to derive and test M 1 and M 2 .
The TOFs are initially calibrated for range-offset and crosstalk to operate through the PDMS. The TOFs are further calibrated with a webcam (C922, Logitech, Lausanne, Switzerland) by implementing Python OpenCV CSRT tracker algorithm [see  the peristaltic actuation protocol. The time-series measurements are compared to determine additional offsets and scaling factors.
The measured RoSEv2.0 conduit displacement data at the center of L 5 , L 6 , and L 7 air pressure chamber of side S 1 are considered as the states of RoSEv2.0. Hence, three discrete states x k,1 , x k,2 , and x k,3 corresponding to layer L 5 , L 6 , and L 7 are defined as RoSEv2.0 states [see Fig. 2(i)]. Each state represents the x-axis conduit displacement at L 5 , L 6 , and L 7 of S 1 , respectively. For modeling simplification, it is assumed that the conduit thickness t remains unchanged throughout its deformation. Additionally, it is also considered that the deformation of the adjacent chambers in each layer is symmetric, and displacement along yand z-axis is negligible [25].

IV. DISCRETE-TIME SINDYC
The original SINDYC algorithm can be extended to discretetime dynamical systems, represented by (1) [28]. There are two primary reasons for implementing the DTSINDYC over the conventional one. First, the calculation of a derivative from noisy data is not required for the DTSINDYC algorithm, which is a significant advantage over its continuous-time counterpart. Second, a discretized model of RoSEv2.0 is necessary for the implementation of MPC.
The data collected from the VPSs or the TOF (x k,1 , x k,2 , and x k,3 , [see Fig. 2 Variables m, n, and l represent the number of samples, state variables (x k,1 , x k,2 , and x k,3 , n = 3), and inputs, respectively, and In this research work, the number of inputs is considered equal to the number of outputs (l = n), where inputs and outputs are digital pressure commands to the valves and TOF or VPS measurements, respectively. The DTSINDYC algorithm implements a library of candidate functions matrix Θ of order m × p, consisting of constant and polynomial (of different order) terms. The library matrix Θ can be defined as where X m 1 ⊗ U defines all the possible product combinations of the components of X and U. Each column of Θ(X m 1 , U) represents the potential candidate for the final right-hand side expression for (1). There is an enormous opportunity for decision in building the library function of nonlinearities. Thus, the system in (1) can be written as where Υ = [υ 0 , υ 1 , . . . , υ n−1 ] is a sparse matrix of coefficients of order p × n such that the sparse vectors {υ i } n−1 i=0 lie in the subspace of R p . The sparsity of υ i selects the active terms from (1). The DTSINDYC model reveals the underlying dynamics of Ro-SEv2.0 in the form of DTDEs. Since only a few of the candidates will appear in the row of g in (1), a sparse regression problem can be defined and Υ can be determined by optimizing it [23]. After the determining Υ, DTSINDYC model of RoSEv2.0 can be written as The MPC optimization is performed over a receding prediction horizon by forward shifting the prediction horizon at every time step N p . At jth time instant, the MPC applies a reference signal x to the cost function minimization (CFM) block (see Fig. 3). The CFM block minimizes a cost function J over N p . The output of the CFM is either used as an input to the valves or DTSINDYC model (M 1 or M 2 , see Table I).
Since the model in MPC is used for prediction; thus, the model output at any instant j is denoted byx j . Between j and j + 1 time step, the single-pole double-throw switch Sw 1 is set to the DTSINDYC model whose output is used to calculate a set of optimal control values u .|x j := {u j+1 , . . . , u j+k , . . . , u j+N c } over a control horizon N c . The behavior of RoSEv2.0 over N p is estimated by its model. To determine the sequence of estimations {x j ,x j+1 , . . . ,x j+k , . . . ,x j+N p −1 }, Sw 2 is initially set to the RoSEv2.0 TOF or VPS forx j and then to the model output for the subsequent predictions. Next, when the CFM algorithm has solved for the best input, S 1 is moved to the valves and first control value from u j+1 is applied to the valves. At next time instant j + 1, the computation is repeated with horizon moved by one time step. The implicit control law for each time step is given by (2). The cost function optimized by the CFM block at each time step can be written as arg min s.t.x j+k = g(x j+k−1 , u j+k−1 ), u lb ≤ u k ≤ u ub . (7) At every iteration (time step) j, CFM minimizes the functional J in (6), which includes a terminal cost atx N p . The leftmost term in (6) penalizes the root-mean-square deviations of the RoSEv2.0 model predicted trajectoryx j+k [given by (7)], from the reference trajectory x (ref) j+k . Control output sequence u j+k and change in control output sequence (u j+k − u j+k−1 ) are also penalized by the middle and the rightmost terms in (6), respectively. The range of u k at every iteration of MPC is constrained by control bounds [u lb and u ub in (7)]. In (6), the weight matrix Q multiplies the deviation of the controlled variable prediction at j + kth instant from the respective reference point. R u and R Δu penalize the value of the control vector and future control moves, respectively. In (6), Q andQ are positive semidefinite, and R Δu and R u are positive definite matrices.
Closed-loop stability for online MPC is an issue if receding horizon predictive control is implemented in such a way that requires online redesign. In MPC, the terminal cost and terminal constraint ensure recursive feasibility and stability of the closed-loop system, but the terminal constraint reduces the region of attraction. In practice, a sufficient prediction horizon size provides necessary condition for the stability and it makes the system stable and feasible in a (large) neighborhood of the origin. On the basis of horizon length, MPC can be categorized into finite and infinite horizon.
Unlike finite-horizon MPC, closed-loop stability is inherent in infinite-horizon constraint less MPC, since at k + 1th instant, no new information enters the optimization problem, so the optimal trajectory determined at k + 1 is same as that of the tail of kth instant [29].
Theorem 1: Suppose that MPC is obtained for a plant, defined by the nonlinear time-invariant equation in (1), by minimizing the cost function in (6). If the plant model remains unchanged and stable throughout the computation of the control law, and settling to the reference point occurs within the prediction horizon, then despite the use of finite horizon, MPC is closed-loop stable.
Proof: The infinite-horizon formulation of (6) can be written as Let u j+k = u c such that u c is the control input required to fixx j+k at x 0 and hence, u c = 0 is required to drive the model output to zero. Therefore, for k > N c , only first term of (8) given by ∞ k=N c x j+k 2 Q remains, which is a matrix geometric mean that converges if the plant is stable. In the vicinity of the origin and zero input, a nonlinear system, such as (1), can be approximated as that of a linear system x k+1 ≈ Ax k . If the linearized system is stable, then all the eigenvalues of A lies inside the unit circle and, hence, In (9), letQ = ∞ i=0 (A T ) i QA i , thenQ satisfies the discrete-time Lyapunov equation [30] given by In (10),Q is a positive semidefinite matrix, if Q is positive semidefinite. Therefore, (8) can be rewritten in a generalized form as When x (ref) = 0, (6) and (11) are similar, where the predictive control problem in (11) is formulated over N c .

VI. MPC DESIGN AND IMPLEMENTATION IN ROSEV2.0
This section ties together Sections III-V in order to design and implement the MPC in RoSEv2.0. In the domain of mathematics, the biological peristaltic waves are typically modeled as per the sinusoidal waveform given by (12) [31], [32] In (12), is minimum radius (mm), z is axial displacement (mm), t is time (s), H(z, t) is conduit's radius as a function of x and t, a is amplitude of the wave (mm), c is velocity of the travelling wave (mm · s −1 ), and w is wavelength (mm). Equation (12) is a modified version of the original equation described by Dirven et al. [27] and Misra and Pandey [31] to incorporate the entire ascending half-cycle (full wavelength). By considering z = z d , (z d is the axial adjacent chamber spacing, which is 15 mm), t d = z d /c continuous time frequency f d = c/w, and t = kT s , discrete-time form of (12) can be written as where discrete frequency ω = 2πf d /F s , sampling frequency The peristaltic waveform profile, given by (13), can be accomplished by controlling time-varying pressure trajectories into the adjacent whorl of chambers simultaneously (see Fig. 4). The pneumatic valves associated with each layer of RoSEv2.0 command the pressure profile in each layer. The role of the MPC is to find out the optimal control law in terms of digital values, which are supplied to the valves to achieve peristaltic wave shapes, defined by the parameters (c and w/2) of (13). Next, three time-shifted versions of (13) for layer L 5 , L 6 , and L 7 of RoSEv2.0 are formulated (x  (14) x (ref) In (14), j,n = H j,n |n = 1, 2, 3}, (14)] and a set of future reference signal values and feedback value ({x j−1,n |n = 1, 2, 3}) from TOF sensors or VPSs located at layer L 5 , L 6 , and L 7 of Ro-SEv2.0 (see Fig. 4). The blocks within the red dashed section of Fig. 3 represent the expanded MPC loop in Fig. 4. By following the minimization procedure discussed in Section V, the MPC loop generates three control values u j+1,1 , u j+1,2 , and u j+3,3 , which are stored in the L 2 , L 3 , and L 4 index of the accumulator. L 1 index is not used since bolus feeder pipes are inserted in RoSEv2.0 conduit and the pipes coincide with L 1 and L 12 conduit locations. Hence, to keep the orientation of the feeder pipes unchanged, layers L 1 and L 12 are not actuated.
In the next iteration (j + 1), a new set of control signals (u j+2,1 , u j+2,2 , and u j+2,3 ) are generated and the earlier set of signals are right-shifted by three in the accumulator. It takes 3q + 1 shifts (iterations) for u j+1,1 , u j+1,2 , and u j+1,3 to reach L 5 , L 6 , and L 7 , respectively, where q is the number of time steps between the peaks of consecutive reference signal given by discrete-time delay q = t d /T s of (13). If c = 20 mm/s, x d = 15 mm, w = 120 mm, and T s = 0.1 s, then t d = 0.75 s, f d = 1/6 s −1 , and q = 7. Here, T s is the time taken by an MPC loop to complete one cycle, which can vary upon conditions, such as reference signal type and model complexity, prediction horizon, and input bounds. With this kind of approach, the MPC has been made robust to adapt for various peristalsis waveform speeds and wavelengths. The only downside to this approach is difficulty in achieving higher wave speeds if T s increases. For best and fast optimization, the MPC was tuned for various prediction horizon lengths, weight matrices, and control input bounds.

A. DTSINDYC Models for Peristaltic Actuation
This section presents the results with the DTSINDYC modeling methodology presented in Section IV. Layers L 5 , L 6 ,  [27]. In this research work, each step of digital value applied to the DACs and valves assembly has generated a corresponding pressure with 1 kPa · step −1 (refer to Section III). The maximum actuation pressure has been restricted to 47 kPa, which is the pressure needed for full RoSEv2.0 conduit occlusion.
A dataset of size 4350 × 3 comprising various amplitude and time-shifted staircase waveform levels [see Fig. 5(a)] and their respective TOF output was collected [orange and blue waveforms in Fig. 5(b) to (d)]. The blue and orange color plots in Fig. 5(b) to (d) also depict the open-loop response of RoSEv2.0, which was collected as a time-series dataset and utilized for training the DTSINDYC model M 1 .
In the first five cycles of Fig. 5(b)-(d), layer L 5 , L 6 , and L 7 have shown displacements up to 10 ± 1.6, 8.2 ± 0.1, and 8.2 ± 0.1 mm, respectively, when actuated with digital DAC values of 38 (or 38 kPa). In the absence of a closed-loop controller, it was challenging to determine the needed digital DAC value for full conduit closure (i.e., displacement of 10 mm) without damaging the robot. The MPC for RoSEv2.0 is implemented to determine the optimal digital DAC values and, hence, pressure values required in the air chambers of each layer to achieve the prescribed peristalsis.
Model M 1 discovered by applying DTSINDYC [33] to the collected dataset is given as ⎡ In ML, more data collection and training does generate an accurate model, but that leads to poor generalization. In such a scenario, the model shows poor performance when presented with new data because too much training data may result in modeling the random noise, referred to as overfitting. It was found that the DTSINDYC model M 1 given by (15) generalizes the entire dataset [black dashed-line plots in Fig. 5(b)-(d)]with a fraction of noisy training data (m = 500 and n = 3), which confirms the robustness of DTSINDYC modeling with low and noisy data availability. Although the model was trained with peristalsis waves of 8-mm amplitude, but it was able to generalize various amplitude peristaltic waves, as shown by the blue color plots in Fig. 5(b)-(d).
Likewise, model M 2 was derived from the data collected from VPSs located at L 5 , L 6 , and L 7 . The recorded output from the VPSs exhibited noisy data in the absence of an output filter. Irrespective of the noise, DTSINDYC robustly predicted the state with great accuracy. The algorithm was able to manage the bias-variance tradeoff adequately to deliver a precise fit to the data. To emphasize the versatility of DTSINDYC, it is important to note that the procedure followed to build and validate M 2 was kept same as M 1 , except that the data for the respective models were collected from different sensors.

B. Stability Analysis
A system with a bounded control input is said to be inputto-state stability (ISS) if each and every state trajectory stays bounded irrespective of the initial state, and if the input grows small, the state trajectory eventually reduces [34]. Rewriting (1) to incorporate the form of difference equations in (15) In (16), g (.) andũ k = [ũ k,1 ,ũ k,2 ,ũ k,3 ] T represent the nonlinear and second and third terms (input terms) of (15), respectively. Equation (16) has an equilibrium point at the origin. The state matrix A determined from the linear approximation of (16) for u k = 0 has eigenvalues at 0.87, 0.82, and 0.91, which proves that the system is asymptotically stable since all the poles are laying inside the unit circle. By definition, a continuous function V on R 3 for (16) becomes ISS-Lyapunov function if it holds the following two properties: , and Property 2: where α 1 (.), α 2 (.), and α 3 (.) ∈ class K ∞ functions and σ(.) ∈ class K function [35].
For another variant of RoSEv2.0, the coefficients of the discovered difference equations given by (15), but a similar kind of analysis can be performed for any set of nonlinear difference equations, provided the equations are asymptotically stable, and inputs are bounded. Intuitively, a soft pneumatic robot is an underactuated and underdamped system, which is always stable under zero actuation. If a unit-step input is applied to the RoSEv2.0, the chamber deformation always remains at the steady-state position, making the robot inherently in ISS.

C. MPC Performance Testing
The results presented in this section are associated with the theory and methodology presented in Sections V and VI. Raspberry Pi 3B+ based host computer (see Fig. 6) is used to run the DTSINDYC-MPC algorithm to calculate the optimal control values required for tracking the desired peristaltic waves.
For implementing the MPC, the Q, R u , and R Δu matrices in (6) were chosen to be positive definite diagonal matrices. By assigning appropriate values to these matrices, the relative importance of the control goals have been considered. In this research work, larger Q values (Q = 5I 3×3 , where I 3×3 represents identity matrix) relative to the R u and R Δu values (R u = 0.5I 3×3 , R Δu = 0.5I 3×3 ) were selected since larger R u and R Δu values result in slower dynamic responses.Q was determined by solving (10).
The control performance of the implemented MPC with model M 1 was estimated in terms of root-mean-square error (RMSE) (mm) and total execution time (s) taken by the MPC algorithm to stop. Both the performance indicators were evaluated by varying the prediction horizon N p from one to eight, keeping the control horizon N c = N p . The evaluation was done prior to the tuning of the control bounds.
To compare the MPC performance of N p > 1 with N p = 1 (single-step ahead prediction), the RMSE for N p = 1 was subtracted from RMSE of all N p . From Fig. 7(a), it can be seen that rmse 1 , rmse 2 , and rmse 3 decrease up to N p = 4, where rmse 1 , rmse 2 , and rmse 3 are the RMSE between the reference displacement signals (x ) and the RoSEv2.0 TOF outputs (x 1 , x 2 , and x 3 ) in response to the optimized control law (u 1 , u 2 , and u 3 ). Since the RoSEv2.0 conduit radius is 10 mm, so a slight increment in the error (> 0.4 mm) will cause a significant impact on the MPC performance. Until N p = 4 in Fig. 7(a), all the errors tend to decrease with an increase of N p . From Fig. 7(b), it can be observed that the execution time drastically goes up as N p increments. Since no significant difference in the MPC performance was found for N p > 4, thus N c = N p = 4 was chosen.
Another reason for choosing N p = 4 is the generation of peristalsis waves in the range of 20 to 60 mm · s −1 , which is the physiological wave speed range in a human esophagus. Since the axial distance between two adjacent RoSEv2.0 chambers (located at L 5 and L 6 ) is z d = 15 mm, thus to achieve c ≤ 60 mm · s −1 , a time difference of t d ≥ 0.25 s is needed between the actuation of the chambers located in L 5 and L 6 . At Np = 4, an iteration time T s = 0.16 s was recorded, thus number of samples N s ≥ t d /T s = 2 were generated. If Np increases, then T s increases, and at T s > 0.25 s, undersampling occurs, and it becomes impossible to generate a wave speed of 60 mm · s −1 .
After selecting a suitable prediction horizon and tuning the control bounds, the performance of MPC with model M 1 was tested on the setup, as shown in Fig. 6. Reference time-series peristalsis waves (x In Fig. 8(a)-(c), it can be seen that the three RoSEv2.0 states (x 1 , x 2 , and x 3 corresponding to L 5 , L 6 , and L 7 , as shown in Fig. 6) tracked the prescribed peristalsis trajectory satisfactorily when actuated with the MPC generated control signals (u 1 , u 2 , and u 3 of Fig. 8(d) and It can also be observed that the accuracy of the tracking is far better while inflating the chambers than deflating, which is not a matter of concern because the inflation of the chambers provides shape to the bolus tail and pushes it forward. Hence,  half reference wavelength was achieved successfully. Peristalsis waves of amplitude and frequency of 10 mm and 140 kHz, respectively, were prescribed to verify the performance of the MPC under higher frequency, which is different from the training data frequency (see Fig. 9). The experimental results from the controller have shown an RMSE of 0.8 mm, which further reinforces the robustness capability of the MPC.
To emphasize the fact that the presented control framework can be extended to any soft-robotic system with any sensor, MPC with model M 2 , controlling the air chamber pressure, was implemented. The MPC with feedback from the VPSs (see Fig. 6) has shown successful tracking of the reference pressure signals (see Fig. 10). In Fig. 10, only state x 1 has been shown. States x 2 and x 3 have shown a similar satisfactorily tracking profile .

D. Application of MPC in Stent Migration Testing
A starch-based fluid thickener (Nutulis, Nutricia, Schiphol, NL) was used to replicate the masticated boluses throughout the experimentation [36]. Boluses were formulated by mixing 1, 2, and 3 scoops with 200 ml of water to cover syrup, custard, and pudding-like consistencies, respectively. The associated concentrations and viscosities are 60, 80 and 100 g · L −1 and 0.45 ±  0.2, 1.2 ± 0.4, and 3.0 ± 1.0 Pa · s −1 , respectively [27]. After the RoSEv2.0 conduit was lubricated with artificial saliva (Aquae Dry Mouth Spray, Hamilton, Reno, NV, USA), the candidate stent was deployed in the RoSEv2.0 conduit by using an intruder sheath and a retrieval thread. The complete stent deployment setup is discussed in [7].
Stent migration experiments with bolus swallow were performed on RoSEv2.0 by using peristaltic actuation protocol with 50 peristalsis cycles. The applied methodology for conducting the experiments has been discussed in [7]. From Fig. 11, it can be observed that with the increase in bolus concentration, the migration slightly increases. The obtained results are in line with the results presented in [7] in which the experiments were conducted in an open-loop manner, and the actual conduit displacement was not tracked. Since stent migration depends on the peristalsis parameters, no significant difference in migration results was found compared to the results presented in [7]. The negligible migration of the stent is due to its exertion of sufficient radial force on the conduit wall, which ensures its proper fixation under peristalsis. But, with the addition of MPC, it can now be assured that prescribed peristalsis has been achieved because in stent implanted RoSEv2.0 experiments, the conduit displacement was controlled by the MPC to ensure occlusion. Hence, implementing MPC has further enhanced the application of stent migration testing in RoSEv2.0.

VIII. CONCLUSION
RoSEv2.0 was developed with embedded TOF sensors to measure its conduit deformation, which has extended its capabilities. Additionally, with the aid of the integrated TOF sensors, the study has first attempted to implement a closed-loop controller in the form of MPC in RoSEv2.0 (see Fig. 4). Since, the methodology presented to design the MPC can be extended to any other soft robotic application provided the SINDYC model can be generated, thus MPC was also implemented to control RoSEv2.0 air chamber pressure. The procedure followed for the implementation remained the same, except for the DTSINDYC models, which were identified on VPS data. Peristalsis waves of speed 20 mm · s −1 , wavelength 75 mm, and amplitudes 5, 7.5, and 10 mm were successfully generated by the MPC (see Fig. 8) for stent testing.
Although the generated model has provided satisfactory performance for peristalsis generation and stent testing, as seen in the results, there is always some model mismatch. Generally, the model mismatch problem in an MPC is treated as a disturbance, which can be addressed by determining the region of attraction (ROA). By introducing a positive invariant terminal set, the ROA can be estimated, where the MPC ensures recursive feasibility and stability. Along with the ROA determination, an online-learning-based MPC in which the model gets updated over time will be designed in future studies to address the elastomer stiffness changes.

ACKNOWLEDGMENT
The code used in this work is made available at: https: //github.com/bhattner143/SINDYc_MPC_RoSE_symmetric _peristaltic.git. The data can be generated using the code.