Edge-Enabled Adaptive Shape Estimation of 3-D Printed Soft Actuators With Gaussian Processes and Unscented Kalman Filters

Soft actuators have the advantages of compliance and adaptability when working with vulnerable objects, but the deformation shape of the soft actuators is difficult to measure or estimate. Soft sensors made of highly flexible and responsive materials are promising new approaches to the shape estimation of soft actuators, but suffer from highly nonlinear, hysteresis, and time-variant properties. A nonlinear and adaptive state observer is essential for shape estimation from soft sensors. Current state estimation methods rely on complex nonlinear data-fitting models, and the robustness of the estimation methods is questionable. This study investigates the soft actuator dynamics and the soft sensor model as a stochastic process characterized by the Gaussian process (GP) model. The unscented Kalman filter is applied to the GP model for more reliable variance adjustment during the sequential state estimation process than conventional methods. In addition, a major limitation of the GP model is its computational complexity during online inference. To improve the real-time performance while guaranteeing accuracy, we introduce an edge server to decrease the onboard computational and memory overhead. The experiments showcase a significant improvement in estimation accuracy and real-time performance compared to baseline methods.


I. INTRODUCTION
S OFT robots made of compliant materials have shown their advantages for adaptability, safety, and novel actuation [1].Nonetheless, the continuum and hysteresis properties of soft robots pose a challenge in analyzing their movement and kinematics compared to the dynamics of conventional rigid robots [2].To address this issue, researchers study and apply continuum mechanics to analyze the kinematics of the soft material, primarily using the visco-hyperelastic model [3], piecewise constant curvature models [4], and finite-element method (FEM) [5].Those kinematic analysis studies establish the soft actuators' plant models and are used to predict the shape deformation under external actuation.A high-precision plant model is essential for the accuracy of robot's actuator control.However, given the gap between the model and the real physical device, an increasing uncertainty will evolve if only the state prediction model is used.Therefore, soft sensors that measure the shape of the soft actuators and have minimal interference to their free movement are emerging and gaining interest within the field of soft robotics [6].
3-D printing has been intensively applied to fabricate soft actuators directly integrated with printed soft sensors [1], which either use commercial conductive filaments [7] or self-developed materials [8] that exhibit deformation-responsive properties.These sensors usually have piezoresistive or capacitive effects, which provides feedback through the resistance or capacitance variations.With the external electrical circuit, sensor feedback can be captured to reflect the deformation of the soft actuator.Such a soft actuator-sensor system establishes a state estimator [9] and enjoys the proprioception property [10], i.e., the system itself has the sensing capacity without the necessity of external sensors, e.g., a camera module with image feedback [11], where its detection accuracy can be largely influenced by the line of sight of the camera in a multiactuator intersection scenario.Thus, the applicability and flexibility of the actuator system are enhanced.
Despite the advantage of proprioception, piezoresistive and capacitive materials cause soft sensors to manifest substantial hysteresis and nonlinearities [12], which confine their usage since basic calibration techniques with regression cannot adequately identify their features.Compensation strategies are desired to resolve this issue.Operator-based techniques identify the hysteresis behavior using the hysteresis operator [13], [14], [15], which models hysteresis as a cumulative memory effect of delayed relay components.These models share a similar operating principle, in which the parameters of highly nonlinear basis functions are fitted through a complex identification process.Furthermore, hysteresis modeling of soft sensors also utilizes data-driven methods.The dynamics of the soft system is treated as a black box model, where the hysteresis is regarded as an inherent property of the sensor [16].
A general limitation of the aforementioned methods is that they essentially construct a deterministic model to characterize the system dynamics, which fails to consider model uncertainty or environmental noise.Moreover, the soft sensor is susceptible to external noise [12], and the dynamics drift in the soft actuator due to polymer aging [17] is likely to be nonnegligible.To overcome these problems, a probability distribution model can be used to characterize the entire soft actuator-sensor system.Considering the dynamics of soft actuators and sensors as a random process allows the state estimation problem to be turned into a Bayesian filtering problem, but only a few studies [18], [19], [20] have been conducted in this area [19], mainly owing to the modeling complexity of nonlinearity and hysteresis with probability models.In addition, during the filtering process, an extended Kalman filter (EKF) linearizes the nonlinear dynamic model, but it is challenging to guarantee system approximation accuracy with the highly nonlinear system.
This study enhances our previous result [20] by characterizing the soft actuator-sensor system with a more accurate representation.Specifically, the Gaussian process (GP) model is employed to identify the nonlinearity and hysteresis in the soft material.The Bayesian nonparametric formulation of GP makes it ideal for modeling the uncertainty and noise in the dynamics of soft sensors and actuators.Utilizing GP regression, researchers can construct the time-serial dependence in soft materials [21], [22], [23].Moreover, owing to its natural Bayesian interpretation and predictive variance, the GP model can also be used to combine the state transition model and the measurement function into the Bayesian recursive estimation [19].By utilizing GP regression, the Bayesian estimation approach can automatically learn the model and noise process from the training data.In light of this, we incorporate the GP models for the soft actuator and sensor into the unscented Kalman filter (UKF) to perform the recursive state estimation.The UKF applies the unscented transform to realize a more accurate approximation than the EKF [9] and exhibits enhanced accuracy for highly nonlinear functions.
The integration framework of GP and UKF was initially proposed by Ko et al. [24].Recently, it has been applied in various fields, such as grip force estimation [25], rock movement prediction [26], and lithium-ion battery state of charge estimation [27].However, its application in soft robotics is still unexplored.Specifically, the hysteresis property induced by soft materials brings high nonlinearity and strong time dependence, which challenges the GP model training and its integration with the UKF.A standard GP-UKF approach cannot accurately estimate the states of soft robots, and the tailored design of the system model for hysteresis characterization of a soft actuator-sensor system remains an open issue.
Another limitation of the GP approach is its computational complexities for model training and online inference, which are cubic and quadratic in the number of training points, respectively [24].The increase in the number of training points also raises a memory burden during online implementation, since the storage capacity of a microcontroller is usually limited.To resolve this issue, we resort to the concept of edge computing, which offloads the computational and memory overhead to the edge server [28].We set up wireless communication and only transmit a small amount of necessary data for executing computationally heavy estimation tasks, by which the communication latency is reduced.With the offloading scheme, the real-time performance of the entire solution is improved without any significant loss of accuracy.
The contributions of our article are summarized as follows.
1) A GP model is formulated to identify the dynamics of soft actuators and soft sensors, which can characterize the hysteresis and uncertainty in the soft material without using parametric models.Consequently, the GP model is more flexible and not limited by the conjectured parametric function.2) We integrate the GP model with the UKF and implement the GP-UKF approach to improve the state estimation accuracy.The state-space model for the soft actuator-sensor system is designed and adapted to characterize the time dependence induced by hysteresis.The developed method is evaluated in real time to test the estimation precision.
The GP-UKF showcases the advantages compared to the pure GP and the previous multihypothesis extended Kalman filter (MH-EKF) approach in estimation accuracy and computational resource efficiency.3) To circumvent the problem of high computational and memory overhead in the GP inference, we offload the online estimation task from the microcontroller to an edge server.The offloading architecture guarantees real-time performance and high estimation accuracy.The edge computing scheme via wireless communication provides more flexibility for future applications of the soft actuator.The rest of this article is organized as follows.Section II first constructs the probabilistic model to represent the hysteresis, with which the dynamics models of soft actuators and soft sensors are identified.Then, the GP-UKF approach is derived and used as a filtering approach to improve the estimation accuracy.We introduce the experimental setup and offloading scheme in Section III.Section IV shows the evaluation and results of our approach.Finally, Section V concludes this article.

A. System Modeling by the GP
This article uses the same type of soft actuator-sensor system in our previous study [20], which consists of a soft actuator and the proprioceptive soft sensor, as illustrated in Fig. 1.The soft sensor is printed onto the soft actuator directly via multimaterial 3-D printing methods.Its electric resistance is deformation responsive owing to the changed distribution of carbon black during deformation.This proprioceptive soft actuator-sensor system constructs a state estimator.
We first elaborate the dynamic models of the soft actuator and sensor, whose state variables are x a,t and x s,t , respectively.The subscript a/s represents actuator/sensor, and subscript t represents the time step.Inspired by the probability box method proposed in [19], we formulate the hysteresis of the soft actuator and sensor by considering their output defined by a probability distribution function.The sensor state x s,t consists of three elements.To express the time dependence, we consider that the sensor output is influenced by values of the bending angle in two consecutive periods of time (i.e., θ t−1 and θ t ).In addition, owing to the memory effect of hysteresis [12], we include the sensor value of resistance in the last time step r t−1 as a state variable, which leads to the sensor model function f s : R 3 → R, and Similarly, the soft actuator f a : R 3 → R can be modeled with where P t is the input pressure.Due to the reality gap, the model functions f s and f a have modeling errors.Therefore, this study applies GP models to describe the probabilistic distribution of functions f s and f a .Due to their similarity, we only show the derivation of GP regression of the soft actuator as an example.The GP model of the soft actuator is characterized by the mean function m a (x a,t ) and covariance function k a (x a,t , x a,t ), i.e., f a (x a,t ) ∼ GP(m a (x a,t ), k a (x a,t , x a,t )).In the formula, where t and t are two time indexes.The mean function is normally initialized as zero with m a (x a,t ) = 0.Moreover, this study selects the squared exponential kernel as the kernel function for the covariance function, with the expression of [29] where λ 2 f is the amplitude hyperparameter denoting signal variance, λ 2 ω is the hyperparameter representing the measurement noise variance, ) is the scaling factor along each dimension of the input space, and δ(•, •) is the Kronecker delta function.The kernel function in (3) follows the form by considering the noisy observations (see, e.g., [30,Sec. 2.2]).The kernel hyperparameters λ λ λ = [λ 2 f , λ 2 ω , Λ] have a significant impact on the shape of regression functions, and they can be optimally defined by log marginal likelihood [29].With a training dataset T a = {X a , θ θ θ} of n training input samples X a = [x a,1 , . . ., x a,n ] and their corresponding outputs θ θ θ = [θ 2 , . . ., θ n+1 ], the optimal values λ λ λ * a can be acquired by maximizing the likelihood of the observed outputs through where K a is the n × n covariance matrix and Note that K a is symmetric and positive semidefinite because of the influence of λ 2 ω .Then, given an arbitrary test state x a,t , the function value f a (x a,t ) is jointly Gaussian with the training inputs X a , with the prior as where ) is an n × 1 vector, and each entry denotes the kernel function value between x a,t and x a,t ∈ X a .0 is a column vector of zeros with the length n + 1, as we consider the prior mean function as zero (i.e., m a (x a,t ) = 0).Based on (6), the distribution for the predicted output θ t+1 is Gaussian, yielding and where μ θ,t+1 and σ θ,t+1 are the mean and variance of θ t+1 , respectively.The sensor model estimation f s through GP regression can be derived similarly.We omit its derivation here for simplicity but will give the necessary formulations in the next section when integrating with the UKF.

B. UKFs With GP Models
Note that when we derive the GP formulation in the previous section, the input state x a,t is regarded as deterministic points without considering uncertainty.Nonetheless, the GP model outputs (i.e., θ t+1 ) are expressed in the stochastic Gaussian Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
form.During the online estimate where the long-term prediction is performed, if we directly incorporate the stochastic state θ t generated by the GP model into the state expression in (2b), the nondeterministic state x a,t will become a non-Gaussian prediction of θ t+1 during the iterations.Such implementation will be problematic in the long run, as already demonstrated in previous studies [19], [31].Kim et al. [19] circumvent the variance and only utilize the mean value from the GP model, by which x a,t can be expressed deterministically.However, the method does not consider the distribution of the input signal and may lead to a growing estimation error.In addition, Deisenroth et al. [31] propose the assumed density filter, where the Gaussian form is preserved through moment matching, and it can capture the true mean and the true covariance of the posterior distribution [9].However, its high computational complexity restricts the real-time implementation.
These problems motivate us to investigate a new method to perform the online estimation with stochastic states in the field of soft robotics.Since the GPs naturally consider the uncertainty and noise in the model, we can improve the state prediction accuracy through the sequential state estimation.Specifically, this study leverages the UKF to retain the Gaussian distribution with high efficiency.Compared to the EKF approach [19], [20], the UKF performs unscented transform, which is a stochastic approximation by a weighted statistical linear regression process.In addition, the UKF is a derivative-free filter that does not require the computation of Jacobians.By contrast, EKF implementation will bring difficulties in our case, since the dynamic functions (i.e., f a and f s ) are not explicitly expressed in the GP approach.Also, compared to other sample-based filters, e.g., particle filters, the UKF adopts the Gaussian distribution assumption, which matches the GP posterior output and improves the computational efficiency in real-time applications.
In this section, we present the unscented transform by referring to the definitions in [9] and [24].Same to the last section, we continue to use the actuator model as the representative.In addition, note that the GP models derived in (1a) and (2a) are regarded as the state prediction function in the Kalman filter, and we shall denote their outputs by rt and θt+1 in the following.
To fit in the UKF scheme, we express a nonlinear actuatorsensor state-space model in the discrete-time domain by where f and g denote the state transition and observation model function, x t is the system state at time step t, and u t is the input signal.z t+1 is the observation value, which differentiates rt since the latter is the ideal estimation.We regard the soft actuator-sensor system as an entirety and define the state by x t = [θ t r t P t−1 ] and the input by u t = P t to reflect the hysteresis influence.The state x t is viewed as a random model, which is represented by the posterior probabilities based on the past measurements and controls.We denote by θ t ∼ N (μ θ,t , σ θ,t ) and r t ∼ N (μ r,t , σ r,t ) the estimated bending angle and sensor value in x t .Thus, the state is only during initialization (i.e., t = 1).
The off-diagonal elements are nonzero after the measurement update in the UKF.Moreover, considering the influence of external disturbances, the input pressure is also stochastic with time-invariant variance σ P .Thus, the system model ( 9) can be expanded with where θt+1 is the estimation value from (2), C = [0 1 0], and z is the sensor measurement noise with z ∼ N (0, σ z ).
The UKF applies deterministic sampling to build a small set of sample points and propagates them through the nonlinear model functions, from which the new Gaussian estimation is formed.Given the posterior distribution from the time step t, we first extract sigma points from this distribution and pass them through the model transition function f .The sigma points are located at the mean and symmetrically along the main axes of the covariance.The dimension of the sigma points depends on the number of random variables.Since the uncertainty is considered in u t , note that herein the dimension of multivariate Gaussian distribution is n =| x t | + | u t |= 4. We denote q t = [x t ; u t ] ∈ R 4 to represent the random variable in unscented transform.There are total 2n + 1 sigma points where μ μ μ q t = μ μ μ t ; P t , Σ q t = diag(Σ t , σ P ) is a block diagonal matrix, (•) i selects the ith column of the matrix, and λ is a scaling parameter that determines how far the sigma points are spread from the mean [32].The sigma points are given as the inputs for the model function f , and the corresponding output is the set Xt+1 = {x t+1 follows Gaussian distribution and t+1 and Σ[i] t+1 are estimated from ( 8) and (10a).Since the set Q t is transferred through the nonlinear function f , the output Xt+1 will not likely be distributed in the Gaussian form.In order to possess the Gaussian distribution property, we extract the mean μ μ μ t+1 and variance Σt+1 of the assumed Gaussian from the mapped sigma points Xt+1 with Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where W m = {w c | i = 0, . .., 2n} are the weights for sigma points when recovering the mean and covariance of the Gaussian distribution (see, e.g., [9,Sec. 3.4.1]).Equations ( 11) and ( 12) establish the essence of unscented transfer that preserves the Gaussian property after the nonlinear function computation.Instead of performing linearization around the mean point as in the EKF, the unscented transform attains the state prediction with Gaussian form without the calculation of Jacobians.
By performing unscented transfer based on ( 11) and ( 12), we obtain the prior distribution of the state at time step (t+1) following Gaussian distribution with xt+1 ∼ N (μ μ μ t+1 , Σt+1 ), where xt+1 = [ θt+1 rt+1 P t ] .Note that θt+1 is different from θt+1 , where θt+1 is obtained from a single point (i.e., μ μ μ t ) projection through f a , while θt+1 is the unscented transform on f when taking the posterior distribution p(x t ) as the input, which enables us to capture the overall uncertainty during the state prediction.
We use the observation model g and measurement z t+1 for the state estimation update in the following step.Note that in (10b), the observation model is linear; thus, we can skip the unscented transform on g as in the conventional UKF, and the predicted observation zt+1 ∼ N (μ z,t+1 , σz,t+1 ) can be expressed directly with Accordingly, the Kalman gain K t+1 ∈ R 3 is calculated by Subsequently, with the measurement z t+1 , the posterior state estimation x t+1 ∼ N (μ μ μ t+1 , Σ t+1 ) can be updated by Thus, shape estimation p(θ t+1 ) can be derived from (15).Note that the off-diagonal elements in Σ t+1 are nonzero during update due to (15b).Algorithm 1 presents the pseudocode for the GP-UKF online implementation.The initialized variance values in Line 1 are estimated based on the observations and experiment, and more details about the numerical parameter initialization are presented in Section IV.Note that lines 2 and 3 can be estimated offline based on the collected training data before the real-time implementation.
As mentioned previously, one main advantage of integrating GP models with the UKF is that the values of estimation uncertainty (i.e., ) can be calculated adaptively by leveraging GP regression, whereas in the conventional UKF, they are manually tuned and often static.With a higher density of training data around the model state, a lower level of state uncertainty variance can be determined, thus achieving the automatic adaptation of the noise variance value.

A. Data Acquisition System
A thermoplastic polyurethane (TPU)-printed soft actuator is fabricated for the validation of the proposed method.Meanwhile, we print the soft sensor onto the actuator directly via multimaterial 3-D printing methods with NinjaTek EEL TM material, which is a mixture of TPU and carbon black [33].
Fig. 2 gives an overview of the experiment setup.The soft actuator deformation is induced by the air pressure from the proportional pressure regulator (Festo VEAB-L-26-D9-Q4-V1-1R1), which is connected to the pressure source of 6 bar at maximum and controlled by an Arduino Due microcontroller board.In addition, to guarantee the accuracy of pressure outputs from the regulator, we calibrate and verify the pressure output with a pressure sensor (Festo SPAU-P6R-W-Q4D-L-PNLK-PNVBA-M12U) before the experiments.We design a resistance meter circuit using a Wheatstone bridge to convert the varied sensor resistances during bending into analog voltage signals, which are then read by the microcontroller board as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the sensor measurement.A compiled program is downloaded to the microcontroller, which performs the estimation and control tasks.We design a proportional-integral (PI) controller with the feedforward part to realize the bending angle control [20].After receiving the control signal calculated by the microcontroller, the voltage amplifier module sends the corresponding electric voltage level signal ranging from 0 to 10 V to the pressure regulator to drive the soft actuator with deformation.
In addition, we monitor the system status and send reference control commands to the board in real time through MATLAB Simulink external mode.As illustrated in Fig. 1(a), we use a web camera to detect the bending angle of the soft actuator by the image feedback module, which is regarded as the ground truth signal for comparison.The image feedback module uses a Logitech C270 HD web camera, which is located on top of the soft actuator and captures real-time pictures at 30 frames/s.The image capture and processing modules are carefully calibrated before conducting the experiments, and the environment of the camera module is static once we finish the calibration, to ensure that the bending angle accuracy captured by the web camera will not be influenced during the training and testing phases.The gray area on the left of Fig. 2 represents the alternative offloading schemes, where the estimation task is offloaded to an edge server.Its implementation is elaborated in the next section.

B. Offloading Scheme
According to our previous study [20], the step response rise time for the same type of soft actuator is around 800 ms.Thus, we define the estimation and control period in this task as 50 ms to ensure that the control frequency is ten times faster.However, the computational complexity poses a challenge to execute the estimation task with the GP approach on the microcontroller in real time.The microcontroller board has limited capacity in both computational speed and memory size, featuring 84-MHz clock speed, 96-kB SRAM, and 512-kB flash memory.Recall that in (8), the calculation of posterior estimation requires the information of K a ∈ R n×n , X a ∈ R 3×n , and θ θ θ ∈ R n , which leads to the quadratic space complexity O(n 2 ) to the number of training data n during the online GP inference.With the growth of training data, the periodic task may not be schedulable for computation or affordable to write in the on-chip memory.We keep the data as the double-precision type both in the storage and memory to guarantee the accuracy of estimation and bending angle control.Thus, when the number of training points is 64, the memory overhead for the GP parameters alone is over 70 kB, posing a big challenge for SRAM capacity.In light of this, we set the maximal number of training points to 64 in this study.
To further resolve this resource deficiency issue, we resort to an edge computing scheme and offload the estimation task from the microcontroller to the edge.The edge server is a desktop PC with Intel(R) Core(TM) i7-8700 CPU at 3.20 GHz, six cores, and 8.0-GB installed memory (RAM).We build the wireless communication between the microcontroller and the edge server through Wi-Fi.Due to the lack of wireless communication capacity on the microcontroller, we connect it with an ESP82-01S board as the Wi-Fi serial transceiver module, which supports Wi-Fi at 2.4 GHz and up to 72.2 Mb/s of data rate.We select the user datagram protocol (UDP) as the protocol due to its simplicity and efficiency.ESP82-01S board transmits UDP packets between the microcontroller and the edge server through Wi-Fi.
The communication path in the edge-enabled approach can be found in Fig. 3.The microcontroller handles sensor reading and pressure control tasks, where z t and u t are sent to the Wi-Fi module through Serial communication.After that, the data are sent to the edge server by the Wi-Fi module, where the GP/GP-UKF are computed.After the edge computation, the posterior estimation p(θ t ) in ( 15) is then sent back to the microcontroller for shape control in the next period, via the Wi-Fi module for transferring the UDP packets to Serial messages.To avoid the buffering of Serial data, the data transmission frequency is set to the same as the control frequency with 20 Hz.Since the estimation and control tasks are accomplished on different nodes, the end-to-end delay in the closed-loop control system can be attributed to four segments, as denoted by t1-t4 in Fig. 3. t1 and t3 represent the communication delays between the microcontroller and the edge server, which contain both the Serial and UDP transmissions.t2 and t4 denote the computation delays from the estimation and control tasks, respectively.To meet the schedulability requirement, the sum of t1-t4 should be less than 50 ms, thus enabling the latest control signal u t to drive the actuator in each period.Based on the signal ranges, we apply fixed point representation and keep two decimal places for each data point, which limits the communication overhead by sending each data point with 2 bytes.

IV. RESULTS AND VALIDATION
We compare the performance of GP-UKF with the following methods for evaluation.
1) GP [30], [22]: The bending angle θ t is predicted by multistep GP regression through (2) without sensor feedback.2) MH-EKF [20]: The actuator is modeled by a second-order state-space function, and the sensor is modeled by several quadratic functions.The Kalman filter is performed through multihypothesis correspondence.3) Bayes filter (BF) [9], [19]: The models are still represented by GP regression.The posterior distribution is calculated by the BF without the unscented transform.
The BF is only used to compare the variance estimation in Section IV-B.In addition, a web camera is applied to provide the ground truth of the estimation by image processing.We use Camera feedback as the true value to calculate the estimation error.
The training data for GP models are collected by driving the actuator with sine waves with positive offsets so that the minimal value is 0 • .The maximal values of the sine waves change from 10 • to 50 • , and the frequency of all waves is 2 rad/s.The hyperparameters of GP models are identified offline through optimizing (4) with the GPML toolbox. 1 We present the identified kernel hyperparameters of GP models for the actuator (λ λ λ * a ) and the sensor (λ λ λ * s ) in Table I.To initialize the parameters in the UKF, we estimate the values of σ P and σ z through experiments.As can be seen in Fig. 4(a), different pressure references P r are set to the pressure regulator, and the readings from the pressure sensor are recorded to calculate the error.The left vertical axis shows the measured pressure P normalized by the pressure reference P r , and the right axis shows the error calculated by | (P −P r ) P r |.The figure shows that the directly regulated pressure by the pressure regulator is higher for low pressures less than 0.5 bars, but the errors are still at low levels less than 3%.For higher pressures, the errors are always less than 1%.The verification result shows the high accuracy of the pressure regulator.For simplicity, we consider σ P as a constant, and assume the average absolute error sampled in Fig. 4(a) as one standard deviation, by which we calculate σ P ≈ 0.001.
The verification of the output resistance value is performed similarly.We measure its stabilized value when driving the actuator with the pressure of a sine wave multiple times.In addition, we shift the time step by offsetting each driving cycle with the start time at the origin and plot the corresponding sensor output in Fig. 4(b), where the average value is also plotted for comparison.The output resistance value in different driving cycles has a standard deviation of 1.3340, and a variance of 1.7796.Thus, we determine the sensor measurement and initial resistance variance with σ z = σ r,0 ≈ 1.8.

A. Comparison With Benchmarks
As shown in Fig. 5, the shape estimation accuracy is assessed by driving the actuator with different positive reference signals.Note that we keep the maximal bending angles of the soft actuator at 40 • to ensure that it is within the domain of the training dataset.In Fig. 5(a) and (b), we compare the estimation performance with sine and triangle waves of different amplitudes, biases, and frequencies.All three approaches can retain high accuracy with 0 • as the origin.However, performance deterioration can be observed in the MH-EKF with the angle biases, which is an open issue identified in [20] due to the model limitation.Compared to the MH-EKF, the GP-class approaches improve the accuracy by adopting the data-driven method with higher model complexity.In Fig. 5(c) and (d), when the bending angle is static, the GP-UKF outperforms the GP with a more noticeable advantage as a result of incorporating sensor feedback.The reason can be attributed to the imperfect coverage of input space in the training data.Since the training data are composed of sine waves, the pressure is continuously varied, where the stationary input state is not included in the training; thus, the model accuracy is influenced.The GP-UKF has a better performance owing to the additional sensor feedback for the estimation adjustment.
Table II lists the estimation error of different types of signals, where the normalized root-mean-square error (NRMSE) is calculated by the range of data from root-mean-square error (RMSE).The square wave has a much higher max error than the other signals because the response time latency incurs a large error with the sudden change in reference, even though this large error is not that noticeable in Fig. 5(c).The GP-UKF  enjoys a consistent higher accuracy among the three approaches, where the NRMSE value is reduced by 30.937% and 57.614% compared to the GP and the MH-EKF on average.

B. Performance With Different Sampling Density
As mentioned in Section III-B, the GP approach involves a much higher computational and memory overhead compared to the traditional methods.Besides the offloading approach for handling the computationally heavy tasks at the edge, another circumvention is to sacrifice the accuracy by reducing the number of training data and saving computation and memory overhead.We compare the tradeoff between the estimation accuracy and the overhead by checking the performance when the number of training data varies.
Fig. 6(a), (b), and (c) compare the estimation performance under three training data sampling densities, i.e., num64, num13, and num7.As mentioned in Section III-B, the maximal number of training datasets is limited to 64 due to memory capacity, which corresponds to the num64 scenario.The numbers of training points in num13 and num7 are 13 and 7, which are five and ten times less, respectively, compared to the num64 scenario.
The reference signals in Fig. 6 are a random signal and a triangle wave.Note that since two models need to be trained (i.e., actuator and sensor), the amount of training data in the GP-UKF is doubled compared to the GP with the same sampling density level.It can be seen that, with less training data, the estimations of both methods show an increasing deviation from the ground truth Camera curve.Table III gives the estimation error comparison of the random reference signal in Fig. 6.Similar to Table II, the GP-UKF performs better than the GP under different densities.In addition, the GP-UKF at num13 maintains a similar level of accuracy compared to the GP at num64.Even though two separate models are required in the GP-UKF, it can achieve the same level of accuracy when the training data density is five times less, where the memory overhead is only 8% compared to its GP counterpart because of the quadratic space complexity O(n 2 ).This implies that the GP-UKF benefits from sensor feedback and is more resource efficient than the pure GP approach.
We compare the variance estimation from different methods at num13 in Fig. 6(d).The GP-UKF has a consistent higher variance compared to the GP, mainly owing to the unscented transfer contributing to more uncertainties during the prior state estimation.It can be seen that, when the estimation error varies in the num13 scenario in Fig. 6(b), the variance changes accordingly in Fig. 6(d), thus realizing the adaptive variance estimation.To further verify that the posterior belief should have a lower level of uncertainty compared to the prior distribution, we design a BF with GP models and exclude the unscented transfer step.Compared to the GP, the BF has a lower variance value.However, both methods show little ability to adjust the variance adaptively compared to the GP-UKF.

C. Hysteresis Loop Analysis
Since the hysteresis phenomenon exists in both the soft sensor and actuator, we compare the estimated hysteresis loop determined by different approaches as a criterion to evaluate the estimation accuracy.As shown in Fig. 7, we focus on the actuator hysteresis estimation by utilizing the sine wave as the reference signal and drawing the relationship between input pressure P t and estimated angle θ t+1 .Besides, the GP-UKFs with two different sampling densities are plotted separately for comparison.It can be seen that the GP-UKF at num64 characterizes the hysteresis loop with the highest precision, while the GP at num64 has a similar trend to that of the GP-UKF at num13, which verifies the results in Table II.The MH-EKF fails to retain the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.shape of hysteresis loops in multiple scenarios, mainly due to its inability to handle the angle bias.
Fig. 8 illustrates the hysteresis loop of the soft sensor.Only the GP-UKF and the MH-EKF are compared since the GP is excluded from sensor feedback.Fig. 8(a) shows the sensor hysteresis loop in different sine waves of amplitudes and biases, where the GP-UKF has a significant advantage over the MH-EKF.Fig. 8(b) plots the temporal variation of the hysteresis loop owing to the sensor drift [12].It can be seen that the GP-UKF achieves an accurate estimation against the time dependence issue.

D. Computational Efficiency Comparison
Fig. 9 analyzes the relationship between the computational overhead and the estimation accuracy.We use response time [34] to represent the estimation task overhead.The response time of the task is evaluated with both onboard and offloading schemes, where the latter enjoys a higher computational capacity but includes the communication delay.With the increase in the number of training datasets, the values of NRMSE gradually decrease, and a higher estimation accuracy is obtained.Meanwhile, the response time of the onboard approach climbs nearly in a quadratic shape.When the number of datasets exceeds 16, the periodic task of 50 ms becomes overdue for the GP-UKF.Due to its memory limitation, the microcontroller fails to handle the task when the

TABLE IV COMPARISON OF THE EXECUTION TIMES FOR EVERY TASK SEGMENT
number of training datasets exceeds 32.Nevertheless, due to the high-computational capacity in the edge-based scheme, the offloading approach remains schedulable with the increase of training datasets.
As indicated by Fig. 3, the end-to-end delay in the closedloop control system is decomposed into four parts.During the experiments, we compare the piecewise execution time for each part and demonstrate them in Table IV, in which t1-t4 are defined in Fig. 3.The data shown in Table IV are collected by performing the GP-UKF approach in the num13 scenario, where the soft actuator is driven by a random signal.We show data by calculating the average and maximal execution times among the task periods in 5 min.The onboard execution times of the control algorithms are measured by Simulink profiling tools, and the execution times on the edge server are estimated by Python Profilers.Compared t1-t4, the most computationally demanding task is offloaded to the edge server, and the microprocessor only needs to handle the PI control and actuation jobs, which significantly reduces the onboard execution time and ensures real-time efficiency.In addition, the wireless data transmission enabled by UDP brings low communication latency with an average time of t1 3.566 ms and sending time of t3 = 1.449 ms.The communication stability is also examined throughout the experiment by recording the number of lost packets, which results in a packet loss rate of 0.137% and guarantees the control accuracy of the offloading approach.The sum of execution times for all task segments in the closed-loop system is less than one task period, which makes sure that the control signal is calculated and actuated based on the latest sensor feedback and will not suffer any deadline misses.

V. CONCLUSION
This study investigated the probabilistic shape estimation problem in a soft actuator-sensor system.The data-driven model, enabled by GP regression, was applied to characterize both the soft actuator and soft sensor dynamics.In addition, to Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
further improve the accuracy of the bending angle estimation, a UKF was designed and collaborated with the GP models.A higher precision estimation system was achieved with the adaptive variance while addressing the hysteresis issue of soft materials.Moreover, the GP approach required a high computation and memory overhead in the online implementation.To enable real-time performance and guarantee the precision of estimation, we circumvented this issue by leveraging the edge computing scheme.The heavy computational estimation task was offloaded to the edge server through wireless communication, and the edge server approach clearly improved the real-time estimation performance through the verification with different types of reference signals.In future work, we will study soft sensor feedback for external force detection.Also, the drift property of the soft sensor will be characterized.

Fig. 1 .
Fig. 1.Experimental setup.(a) Test rig and 3-D printed soft actuator.(b) 3-D printed soft actuator with the integrated soft sensor.The soft actuator is deformed with the bending angle θ t .

Fig. 4 .
Fig. 4. Variance verification of measurement for input pressure and output resistance.(a) Pressure input verification.(b) Output resistance verification.

TABLE II ESTIMATION
ERROR COMPARISON WITH DIFFERENT TYPES OF SIGNAL