Learning 3D shape proprioception for continuum soft robots with multiple magnetic sensors

Sensing the shape of continuum soft robots without obstructing their movements and modifying their natural softness requires innovative solutions. This letter proposes to use magnetic sensors fully integrated into the robot to achieve proprioception. Magnetic sensors are compact, sensitive, and easy to integrate into a soft robot. We also propose a neural architecture to make sense of the highly nonlinear relationship between the perceived intensity of the magnetic field and the shape of the robot. By injecting a priori knowledge from the kinematic model, we obtain an eﬀective yet data-eﬃcient learning strategy. We first demonstrate in simulation the value of this kinematic prior by investigating the proprioception behavior when varying the sensor configuration, which does not require us to re-train the neural network. We validate our approach in experiments involving one soft segment containing a cylindrical magnet and three magnetoresistive sensors. During the experiments, we achieve mean relative errors of 4.5%.


Introduction
The past decade has seen an explosion of novel continuum soft robotic platforms. [1][2][3] Inspired by invertebrate animals, these robots are almost entirely composed of soft deformable materials. 4 This makes them robust and safe, but at the same time, it renders their modeling, 5 control, 6 and shape sensing 7 substantially more complex than for their rigid counterparts. Shape sensing is especially intricate because it is both a technological and algorithmic challenge. Sensors must not obstruct the natural behavior or reduce the compliance of soft robots. At the same time, non-collocated and nonlinear sensors require algorithms for the measurements to be interpreted and connected to a description of the robot's shape.
Several sensing modalities have been considered to implement shape sensing, such as resistive, 8,9 capacitive, 10,11 optical, 12 and visual. 13 Magnetic sensors [14][15][16][17][18] are a promising solution as they are compact, highly sensitive, and can be easily integrated into existing soft robot designs. Thus, they can provide reliable and fully proprioceptive measurements at the cost of a minimal decrease in the robot's softness. Among the above-cited papers, the only work leveraging external magnetic fields is by Song et al., 14 where coils placed at a distance generate the magnetic field. Along this thought, an obvious choice would be to measure the earth's magnetic field for proprioception purposes. However, it would only allow estimating two rotational components, and any translational effects, such as an elongation of the continuum robot, could not be captured. Alternatively, some papers [15][16][17] use co-axial pairs of magnets and sensors embedded in the robot to estimate planar deformations. Recently, Mitchell et al. 18 have shown that a similar strategy can sense 3D deformations. Such simple arrangements greatly simplify the analysis, allowing for connecting sensor readings to the robot's shape through direct interpolation. Nevertheless, relying on isolated pairs of sensors and magnets also strongly limits the density and the amount of information gathered through this method.
This paper proposes to use permanent ring magnets and multiple magnetic sensors for shape sensing of continuum soft robots. Importantly, we remove the requirement of magnets and sensors being placed in co-axial pairs. We have been inspired by recent work leveraging deep learning to interpret various types of non-magnetic sensor data for proprioception purposes. [19][20][21][22] However, learning end-to-end mappings from sensors to configurations has three significant drawbacks: (i) it is data-intensive, (ii) it calls for recurrent architectures to encourage temporal consistency for the robot's configuration estimates, (iii) it requires re-training when changing the kinematic model of the robot. We propose a neural architecture that circumvents all three issues (see Fig. 1). First, we train shallow neural networks to predict the measurements of the magnetic sensors from a parameterization describing their relative pose with respect to the magnets. We then optimize the configuration estimate -and thus the sensor positions -to minimize the error between the predicted and actual sensor measurements. This way, we introduce a priori information on the modes of deformation of a continuum soft robot, effectively removing the kinematics from the black box. The presented strategy enables us to re-arrange and remove redundant sensors during inference without requiring us to re-train the neural network on the adjusted sensor configuration.
To summarize, this paper contributes to the state of the art in soft robot sensing with: A proprioceptive sensing modality relying on multiple magnetic sensors in conjunction with a neural network-based architecture that learns to estimate the full 3D shape of the robot from the sensor readings.
Injection of kinematic priors through a description to spatially relate the poses of a sensor to the magnets. This proposed parametrization serves as an input to a neural network predicting the sensor measurements.
Experimental verification of the approach for a onesegment robot with three integrated magnetoresistive sensors and proprioception of 3D curvature.
We open-source a Python / PyTorch implementation of the proposed algorithm and the corresponding datasets on GitHub. §

Proprioception with magnetic sensors
This section introduces a methodology to achieve proprioception for soft robots with magnetic sensors. We consider a continuum robot with the shape of its backbone described by the configuration variables q 2 R nq . In the commonly used piecewise constant curvature (PCC) kinematic state parametrization, 23 the continuum robot is assumed to consist of n b segments with each segment exhibiting constant curvature and elongation along its length. Therefore, the configuration of the soft robot can be described with q 2 R 3n b . Please refer to Appendix A.1 for more details. We indeed use PCC for most of our simulations and experiments. Note, however, that the proposed proprioception algorithm applies to any finite-dimensional kinematic description of a soft robot. 5 In fact, we also specifically consider a robot with affine curvature 24,25 with its shape described by the configuration q 2 R 4 . We document this alternative kinematic model in Appendix A.2. The proprioception methodology described in the remaining section is agnostic to the chosen kinematic model.

Proposed method at a glance
We integrate n m axially symmetric magnets and n s magnetic sensors into the robot. The magnets must be installed along the center line of the segment while the sensors can be arbitrarily placed. Fig. 1 concisely describes the proposed shape-sensing strategy with a pictorial example of a soft robot consisting of three segments and equipped with five magnetic sensors and three cylindrical magnets.
The goal of the algorithmic part (center of the figure) is to regress the robot shape (left) by estimating the configuration q of all segments (right), starting from the measurements of the magnetic sensors u (e.g. usually the magnetic flux density). We achieve this by training a sensor measurement predictor and subsequently optimizing the configuration estimate q for the prediction û to match the actual sensor measurement u. Instead of predicting the sensor measurements end-to-end, we decouple the kinematics by explicitly describing the kinematic relationshipx j = f x,j (q) between sensor j and each magnet. Then, we train a neural network f p j to predict the measurement û j based on the inputx j . To achieve proprioception, we jointly optimize q for all sensor measurement predictions. Fig. 1 Proprioception for continuum soft robots with magnetic sensors: an initial configuration estimate q is employed to derive the kinematic relationshipx between a sensor and each magnet. Subsequently, this kinematic description is used to predict the sensor measurementsû 2 R ns with a neural network trained in advance. A mean squared error (MSE) metric evaluates the accuracy of the predictions û compared to the actual measurements u. Finally, we optimize the configuration estimate q to achieve proprioception by minimizing the sensor measurement prediction loss.

Magnet sensor kinematics
This subsection derives the kinematic relationship x j ¼ f x j ðqÞ 2 R 1þ4nm between the jth sensor and all n m magnets as we hypothesize that we can estimate the sensor measurement u j solely based on (a) the angle l j 2 R 1 to the earth's magnetic field, (b) the distance d k j between the jth sensor and kth magnet, (c) the angle a k j between the sensor measurement direction and the cylindrical axis of the magnet, (d) the angle b k j between the sensor measurement direction and the vector from the magnet to the sensor, and (e) the angle y k j between the cylindrical axis of the magnet and the vector from the magnet to the sensor. Accordingly, x j is defined as with x k j 2 R 4 the kinematic relationship between the jth sensor and the kth magnet: We visualize the parameters incorporated in x k j in Fig. 2(a). In the following, we present the derivation of all components of x k j . Please note that all kinematic frames used in the following paragraphs are visualized in Fig. 2 We define that the kth magnet is integrated into the ith segment. Now, we first derive a transformation matrix T m k iÀ1 from the base frame {S iÀ1 } to the magnet frame {S m k }. This can be achieved by evaluating the chosen kinematic model, two of which we report in the Appendix A, at the segment coordinate v ¼ d m k L 0;i . This means that the cylindrical magnet is integrated at a distance, which is measured along the backbone, of d m k from the base of the segment. Subsequently, we describe the pose of the jth sensor with respect to the base of the ith segment. Denote with d s j ,r the radial distance of the sensor from the center line, with j j the azimuth angle of the sensor in the cylindrical plane, and with d s j ,a the axial distance along the center line from the base of the ith segment. We derive the transformation T s j ;r 0 iÀ1 to the center of the cylindrical plane of the sensor analogously as for the magnets by evaluating the forward kinematics at v ¼ d s j ;a L 0;i . This is followed by applying the radial offset d s j ,r in the cylindrical plane of the sensor Optionally, a static rotation offset can be applied to R s j iÀ1 such that local z-axis o s j corresponds to the sensor measurement direction.
Knowing the transformation matrices from the base frame of the respective segment to the sensor and magnet frames, we express them in the inertial frame {S 0 } by multiplying with the kinematic chain T iÀ1 Next, we need to express the sensor measurement direction o s j and the cylindrical axis of the magnet o m k in the inertial frame As the sensor measures contributions of the earth's magnetic field, we need to state the angle l j between o s j and the earth's magnetic field unit vector n e . Similarly, we investigate the angle a k j between o s j and the cylindrical axis of the magnet o m k cos(l j ) = n e Áo s j , cos(a k j ) = o m k o s j .
We define the translation and distance between the magnet and the sensor in the frame {S 0 } as: Building on the derivation in (5), we compute the angles b k j and y k j using the dot product rule Lastly, the kinematic descriptions for all sensors are vertically stacked as We will in the following refer to the mapping f x ðqÞ:q 2 R nq ! x 2 R nsþ4nmns as the magnet sensor kinematics.

Data-driven sensor measurement model
We use a data-driven approach to learn the forward sensor model û = f p (x j ) for each sensor using a neural network parameterized with p. We note that the same neural network weights can be shared for all sensors, but oftentimes performance can be improved by training a specialized model with weights p j for each sensor. During training on a dataset of length n t , we minimize the mean squared error (MSE) error between the predicted sensor measurement û j and the actual sensor measurement u j : where t denotes the current time index. Note that whenever we omit the time index in our notation, we always refer to the current time t. Finally, to simplify the notation, we combine each sensor measurement prediction u j 2 R into an array u 2 R ns and stack the neural networks as f P (x):xu. We discuss the choice of the specific network architecture later in Section 3.

Proprioception by optimizing configuration estimate
Now, that we are able to predict the sensor measurement û using the composition of the kinematics f x (q) and the neural networks f P (x), we need to optimize the configuration estimate q for the predictions û to match the actual sensor measurements u as closely as possible. We capture the error between the predicted sensor measurements û and the actual sensor measurements u by the MSE loss function we strive to minimize Accordingly, the optimal configuration estimate q can be found withq We optimize the cost function (9) through iterative gradient descent, as detailed in Fig. 3. The gradient descent is initialized with the best estimate of the previous time-step q 0 (t) = q(t À 1). As common in literature, we optimize the state belief q l with a step size g and the momentum m using the Jacobian of the loss We can use the chain rule to derive an analytical expression for the gradient of the loss incorporating the gradient of the magnet sensor kinematics q q fx(q) and the gradient of the neural network qxf P (x): After executing the gradient descent for n it iterations, we evaluate which iteration l* had the lowest loss L u and accordingly select qðtÞ ¼q l Ã as the best configuration estimate of time-step t.

Piecewise constant curvature simulations
We evaluate the proposed methodology for estimating the PCC kinematic configuration q 2 R 3n b of soft continuum robots thoroughly in simulations. The PCC model allows for bending and elongation of each segment in 3D space. Please refer to Appendix A.1 for more details. We vary the number of robot segments n b , remove and add sensors (i.e. change n s ), modify the arrangement of sensors, and the direction of the earth's magnetic field n e . To motivate some of the unique advantages of our method, we use the same learned neural network weights for all these trials.

Simulation setup
In our simulations, we consider a robot consisting of one, two, or three segments. All cylindrical segments have an unextended length of L 0,i = 110 mm and a radius of d i = 22 mm. Each segment has a ring magnet of outer diameter 12 mm, inner diameter 6 mm and thickness 6 mm integrated along the backbone at a distance of d mk = 55 mm from the base of the segment. In the nominal case, three sensors are placed symmetrically in the tip plane of each segment at a radius of 13 mm from the center with the sensor measurement direction pointing along the local, negative z-axis of the tip plane {S i }. We also consider alternative placements of the sensors, which we further detail in Section 3.5. We build on Magpylib 26 to simulate the magnetic field behavior. We model the magnets as cylindrical neodymium grade N50 rings with a magnetization of 1450 mT in the local z-direction. After simulating the magnetic field, we rotate the B-field into the local reference frame of each sensor and take the local z-component of the magnetic flux density as the sensor measurement.

Prediction network
The training set consists of 120 000 random configurations sampled from uniform distribution D x;i $ UðÀ20:7 mm; 20:7 mmÞ, D y;i $ UðÀ20:7 mm; 20:7 mmÞ and dL i $ Uð0; 5:5 mmÞ, where the upper bound represents a bending of the tip of 541 with respect to the Fig. 3 A block-diagram of the gradient descent as an iterative update loop for the configuration belief q(t). The gradient descent is initialized with the optimized solution q 0 = q(t À 1) from the previous time-step. Note that during the iterative loop, q l+1 updates q l in the memory block.
base of a segment and an elongation of 5%. We also randomize the placement of the sensors in the training set. While in the nominal case, the first of the symmetrically placed sensors is placed on the local x-axis, we randomly sample an offset angle 2p n s i for each training sample, where n s i = 3 is the number of sensors per segment. Finally, we also randomly sample the radial displacement of the sensors from the center with d s;r $ Uð8:7; 17:3Þ mm and consider a tilting of the sensors (i.e. a rotation around the tangential axis) with c s $ UðÀ20 ; 20 Þ. Before training, we randomly split off 30% of the training set for validation purposes.
We conducted a selection study involving hyperparameters and feed-forward neural network architectures (number of layers, nodes, and nonlinear activation layer types) on the validation set. In particular, we aimed to generate a smooth loss landscape to improve the gradient descent convergence leading us to employ a stochastic gradient descent (SGD) optimizer in conjunction with the stochastic weight averaging (SWA) 27 strategy. The neural network itself has 18 layers in total and contains, after an initial 1D batch norm layer, four blocks and is concluded with a fully-connected layer at the end. Each block consists of a dropout with a probability of 1%, a linear layer, a ReLU, and a 1D batch norm layer. The hidden state is first increased to 96 nodes, then to 256 nodes, and finally reduced again to 64 and 24 nodes. We minimize a mean squared error (MSE) loss of the neural network prediction û j (t) for 250 epochs with batch size 650 while setting an initial learning rate of 0.18 for the cosine annealing learning rate scheduler. 28 The SWA 27 strategy is started after 125 epochs. We train the neural network such that all sensors in the ith segment share the same weights p i . When the training is finished, we select the model from the epoch with the lowest validation loss and save it for later testing.

Optimization
We optimize the configuration variablesq i ¼ ðD x;i ;D y;i ;dL i Þ T for each segment to minimize the sensor measurement prediction loss as defined in (9). The optimization strategy solely relies on gradient descent and uses the best solution from the last time step q*(t À 1) as an initialization q 0 (t). For the first time step of the trajectory, we initialize with the ground truth. We use a step size g = 3.5 Â 10 À4 , 3 Â 10 À3 , and 2 Â 10 À3 during gradient descent for a one-segment, two-segment and threesegment robot respectively. For all robots with more than one segment, the step size is reduced by a factor of ten when optimizing the elongation. The momentum m is 0.3 for all trials and we perform 20 gradient descent iterations for each time step.

Evaluation
We evaluate the performance of our method at estimating the configuration of the segment by computing a relative root mean-squared error (RMSE) metric with respect to the ground-truth configuration q(t) for each configuration variable separately e q? ¼ where n t is the number of discrete time-steps and D q? is the dynamic range of each configuration variable between the maximum and minimum value in the test set. All simulations are evaluated on a full lemniscate trajectory, which is very similar to what is plotted based on experimental data in Fig. 9(d). The maximum bending angle of this trajectory is 45 and the elongation of the segment follows a cosine wave with a minimum and maximum at 1.25% and 3.75% respectively.

Results
We present all simulation results in Table 1. In the first section, robots with one to three segments, which all exhibit a nominal sensor placement, are considered separately. Neural networks are trained separately for each of these robots, as the input dimension needs to be adjusted to the number of magnets. The results show that the method works well for robots with 1-3 segments. It can be observed that the estimation error is usually lower for the distal segment(s). Next, the number of sensors n s is varied for a three-segment robot. We always apply a symmetrical placement of the sensors in the tip plane of each segment. In the first trial, two sensors are mounted in the tip plane of each segment opposite to each other (6 sensors in total). While the bending along D x,i and the elongation of the segments dL i can still be estimated, the setup does not contain sufficient information to accurately determine the bending into D y,i . While the nominal case of nine sensors in total already achieves relative RMSEs in the range of 1.6% to 6%, the proprioception performance can be slightly improved by adding more sensors. We emphasize that the neural networks are not retrained when adding or removing sensors. In Fig. 4(b) and 5, we plot the proprioceptive performance of a three-segment robot with four sensors per segment. Then, we simulated a failure of the 4th, 8th, and 12th sensor at 5s by removing these sensor measurements from (11) (i.e. the gradient descent). Our method is able to adapt without re-training and leverage the nominal redundancy of sensor measurements, as three sensors per segment are sufficient for shape estimation.
Then, two adjusted sensor placements are investigated, while keeping the neural network weights constant. First, the sensors are tilted from the nominal case of pointing along the local z-axis by c s = 101 towards the inside. In a separate simulation, the sensors are moved radially from nominally 13 mm to 16 mm. As the results show, the configuration of all three segments can still be estimated accurately with a mean error of 3.3%.
Finally, an earth's magnetic field of magnitude 0.065 mT is added. A separate neural network is trained on a training set with randomly sampled magnetic field vector directions n e . As the last section of able to adapt to any earth's magnetic field direction by leveraging the l j input parameter.

Affine curvature simulations
To demonstrate the efficacy of the proposed method for higherorder kinematic models than PCC, we also conduct simulations of an affine curvature soft robot. The affine curvature kinematic parametrization 24 has been shown capable of representing the shape of soft tentacles 25,29 and provides a continuous function k = k 0 + k 1 v to describe the curvature of the soft robot, where k 0 , k 1 are two configuration variables and v A [0,1] is the backbone coordinate. We allow for movement in 3D space by also specifying an azimuth bending angle f and the elongation dL. Please refer to Appendix A.2 for more implementation details about the affine curvature model.

Simulation setup
We use the same simulation setup as described in Section 3.1. Therefore, we report in the following only the implemented modifications to simulate an affine curvature soft robot in Magpylib. 26 Namely, we consider one affine curvature segment Sensor measurement predictions (top) and configuration estimates (bottom) for a three-segment robot with twelve sensors nominally. We simulate a sensor failure of the 4th, 8th, and 12th sensors after 5 s of trajectory time by removing the measurements of these sensors from the gradient descent. This can be compensated automatically by the redundancy of the sensor configuration. We plot the groundtruth values u and q in solid, the estimate û and q as a mean over three random seeds with dashed lines and the standard deviation as an error band. Table 1 Simulation results: first, we report the absolute root mean-squared error (RMSE) e u of sensor measurement predictions on the test set averaged across all sensors on the robot. Next, we state the relative RMSE [%] of each robot configuration estimate. All models are trained for each segment separately on a trajectory with randomly sampled configurations and sensor kinematic parameters and evaluated on a lemniscate trajectory. The first section applies our methodology to robots consisting of a different number of segments n b with three sensors attached to the tip of each segment. The number of sensors is varied in the second section for a three-segment robot with all sensors placed symmetrically. The third set of trials then investigates how robust the method is to changing the kinematic parameters of the sensors, such as the tilting angle of the sensors c s and the radial distance of the sensors d s,r . Finally, we apply the earth's magnetic field along different cardinal directions in the inertial frame. The RMSE of the configuration estimates is normalized with the range of the dataset for each configuration variable as stated in (13). We report the error as mean AE stdev and compute the statistics over three different random seeds. The random seed determines at the start of the training the initialization of the neural network weights

Prediction network and optimization
We simulate 120 000 random configurations of the affine curvature robot to generate the training set. For this purpose, we sample the configuration variables from uniform distributions: the affine curvature parameters k 0 2 Uð0; 0:942 rad m À1 Þ, k 1 2 Uð0; 3:770 rad m À2 Þ, the azimuth angle f 2 Uð0; 2p radÞ; and the elongation dL 2 Uð0; 6:6 mmÞ. Before training, we randomly split off 30% of the training set for validation purposes. We train a specialized neural network f p j (x j ) for each sensor and use the same neural network architecture as in Section 3.2 with the exception of the addition of a final layer y(x) = sign(x)e |x| . The training runs for in total 250 epochs and uses the SWA 27 strategy with a learning rate of 0.01. All other training hyperparameters are the same as in Section 3.2.
The four configuration variables are optimized to minimize the loss between the predictions and simulated measurements of the nine sensors as defined in eqn (9). For this optimization procedure, we employ gradient descent running at 40Hz with step sizes of g k 0 = 1, g k 1 = 5, g f = 1, and g dL = 2 Â 10 À4 . The momentum is set to m = 0.3 and 20 iterations are performed at each time step.

Evaluation
We evaluate the trained model on a flower trajectory of duration 10 s and sample rate 40 Hz. The evaluation trajectory has the following characteristics: k 0 is actuated by a sinusoidal wave of frequency 0.3Hz in the range 0; p 4 rad m À1 h i . Similarly, k 1 is also varied through a sinusoidal function of the same frequency and has a dynamic range of [0.1p,p] rad m À2 . The azimuth angle f is linearly scaled from 0 rad to 2p rad over the duration of the trajectory. Finally, dL follows a sinusoidal sequence of frequency 0.1 Hz in the range of [0,5.5] mm. We use the same evaluation metrics as first introduced in Section 3.4. We report the error as mean AE stdev and compute the statistics over three different random seeds. The random Fig. 5 Sequence of stills for simulated sensor failure as shown in Fig. 4(b): the number of sensors per segment, which are rendered in blue, is reduced from four to three at t = 5 s. We visualize the ground-truth shape of the soft robot with full opacity and the estimated configuration with slight transparency. The magnets are rendered in green and the three segments are visualized in a color sequence from black to violet.

Fig. 6
Sensor measurement predictions (top) and configuration estimates (bottom) for an affine curvature segment with nine magnetic sensors. We plot the ground-truth values u and q in solid, the estimate û and q as a mean over three random seeds with dashed lines and the standard deviation as an error band. The random seed determines at the start of the training the initialization of the neural network weights. Fig. 7 Sequence of stills for a simulated affine curvature robot. We visualize the ground-truth shape of the soft robot with full opacity and the estimated configuration with slight transparency. The two magnets are rendered in green and the nine sensors in blue.
seed determines at the start of the training the initialization of the neural network weights.

Results
The trained neural networks achieve an RMSE error for predicting the magnetic sensor measurements of 0.025 AE 0.002 mT on the test set. When we run inference (see Fig. 6) on the flower trajectory, the configuration variables can be estimated with an absolute RMSE of e k 0 = 0.042 AE 0.005 rad m À1 , e k 1 = 0.11 AE 0.04 rad m À2 , e f = 0.08 AE 0.02 rad, and e dL = 0.001 AE 0.001 mm. We state the relative RMSE errors for the configuration estimates as e k 0 = 5.4 AE 0.7%, e k 1 = 4.0 AE 1.4%, e f = 1.3 AE 0.4%, and e dL = 3.8 AE 1.6%. In Fig. 7, we render at six different points along the trajectory the robot's shape according to the ground truth and estimated configurations respectively. The sequence qualitatively shows that our proposed method is able to estimate the affine curvature robot's shape very accurately.

Experiments
We verify the performance of our proposed proprioception method in experiments involving one soft robot segment with three magnetoresistive sensors attached to the tip. We aim to estimate the CC configurationq ¼ ðD x;1 ; D y;1 Þ T 2 R 2 from the measured sensor values uðtÞ 2 R 3 . We let the robot follow a diverse set of trajectories and evaluate the proprioception performance. After measuring the ground-truth pose of the tip of the segment with a motion capture system, we perform inverse kinematics with the closed-form solution reported by Della Santina et al. 30 and quantitatively compare the proprioceptive configuration estimate of the segment with the groundtruth configuration.

Robot design
We use a cylindrical, pneumatically-actuated soft robotic silicon segment of length L 0 = 110 mm and radius d 1 = 22 mm consisting of three independently inflatable cavities evenly spaced in the radial direction from the center line. 31 The proprioceptive sensing system is achieved by embedding one ring magnet in the backbone at a distance d m 0 = 55 mm from the base of the segment and three symmetrically-placed Magnetoresistive Sensors (MRSs) at the tip of the segment as visualized in Fig. 8. Although we use MRSs in our experimental setup for their high sensitivity, 32 this is not a strict condition and other sensor types measuring the magnetic field such as Hall-effect sensors can be combined with the methodology proposed in this paper too. For casting the silicone segment, we use a 3D-printed mold with a holder for the magnet which keeps it in place inside the segment. The magnet used is a neodymium ring of grade N50 with a thickness and inner diameter of 6 mm each, and an outer diameter of 12 mm. The MRSs of type Honeywell HMC1021Z are integrated into a printed circuit board (PCB) and output a voltage difference of 50 mV mT À1 . The three sensors are equally spaced at 1201 from each other and are placed at a radial distance of d s,r = 13 mm, and at a longitudinal distance of d s,a = 116 mm from the base in a straight configuration. For each sensor, we implemented a Set/Reset and an amplification circuit on the PCB. The Set/ Reset circuit is used for calibration of the sensor by re-aligning the magnetic domains. After amplification of the sensor output by a factor of 100, the output of the sensors is processed with a Texas Instruments ADS1115 module resulting in a digital signal of 16bit resolution. All sensor measurements u are in the range [0 mV,2048 mV], which corresponds to magnetic flux densities of [0 mT,41 mT].

Experimental setup
We conducted our experiments in a lab environment with the base of the soft robot segment mounted in a tip-down configuration to a cubical cage as shown in Fig. 8. A 3D-printed ring with four reflective markers is mounted on the tip of the segment. Eight motion capture cameras are attached to the cage tracking at 40 Hz the 3D pose of the ring. We transform the pose measurements of the tip to the base frame of the robot and compute the closed-form inverse kinematics 30 to receive a ground-truth configuration estimate qðtÞ 2 R 2 . Each of the three pneumatic chambers of the segment is connected via tubing to a separate valve of a proportional pressure regulator operated at 100 Hz. We read out the analog signals of the magnetoresistive sensors with an Arduino Uno at 40 Hz and save them for later offline processing. We temporally align the motion capture and the magnetic sensor data by detecting the initial extension of the robot with a suitable threshold. The sensor noise is determined for both an unelongated straight configuration, and during fully inflated bending. Here the standard deviations of the white noise are 0.24 mV and 3.55 mV, which normalizes to 0.03% and 2% of the dynamic range respectively. Furthermore, we identify the earth's magnetic field direction in the base from as n e = (À0.311,À0.234,0.921) T using a compass and the world magnetic model (WMM). 33

Pneumatic actuation and trajectories
We consider, as visualized in Fig. 9, six continuous actuation sequences in this paper: random configuration way-points which are connected through linear interpolation (T0), planar side bending (T1), the tip following a half lemniscate (T2) and full lemniscate (T3), a spiral with constant linear velocity 34 (T4) and finally a flower-shape (T5). We define our trajectories as wrenches t xyz = [t x t y ] T on the tip of the segment in Cartesian space, where t x and t y cause bending around the local x-and yaxis of the tip respectively. The pressures we command from the pressure regulator are given by inversely evaluating the force produced at the center of pressure at the tip of the segment for each chamber for a given chamber pressure. 35 All actuation sequences are preceded by first applying an offset pressure of 225 mBar in all chambers, which causes a nearconstant elongation of the segment. The peak pressure, which causes maximum bending, is set for all trajectories to 450 mBar. The 1D bending (T1), half lemniscate (T2) and full lemniscate (T3) are all executed periodically with a period of 5 s, 10 s, and 10 s respectively. Trajectories T0 and T4 are characterized by a constant velocity in torque-space of 0.025 N m s À1 and 0.0125 kN m s À1 respectively. The flower trajectory T5 can be described as periodic 1D bending with a linearly changing azimuth angle. It exhibits an angular velocity of 0.0126 rad s À1 and a period of 10 s for the bending, which results in 50 bending cycles per circumnavigation. While the random configuration set points of T0 are recorded for 200 s, T1, T2 and T3 have total a duration of 90s each, and the spiral T4 and flower T5 last for 120 s and 1500 s respectively. We split off the final 20% of all datasets as a test set.

Prediction network and optimization
We use the same neural network architecture and training procedure as in Section 3.2, but with an adjusted initial learning rate of 5 Â 10 À5 and train the model separately for each sensor on the 1200 s long T5/flower trajectory, which results in 48000 training samples for each neural network.
We optimize the configuration variables D x and D y for the one segment to minimize the sensor measurement prediction loss as defined in (9). The optimization strategy solely relies on gradient descent running at 40 Hz with a step size g = 1.5 Â 10 À8 and momentum m = 0.2. We visualize a sample loss landscape in Fig. 10(a).

Results
First, we quantify the performance of the neural network predicting the sensor measurements û(t) for a known, ground-truth configuration q(t), which we report in Table 2. We observe, that the relative RMSE lies between 0.6% and 4.6% of the range of the respective datasets with a mean of 3.1%. As expected, the predictions are generally the most accurate when evaluated on a trajectory of the same type as the neural network was trained on (T5). Next, we analyze the proprioception performance on the same trajectories. We report relative RMSEs between 1.9% and 13.6% for estimating the bending of the robot. Additionally, we visualize the configuration estimates for two  With the hue, we visualize the RMSE of the sensor measurement prediction û for a given configuration q = (D x ,D y ) T . Additionally, we denote the initial configuration estimate with q 0 , the trajectory of the gradient descent with q l , the optimal configuration with q, and the ground truth with q. Panel (b) top: Proprioception on the test set of T2 using a model trained on T5. Panel (b) bottom: Configuration estimates for a model trained and evaluated on separated parts of trajectory 5. We plot the ground-truth configuration q in solid, the estimate q as a mean over three random seeds with dashed lines, and the standard deviation as an error band. The bottom right plot zooms onto a selected part of T5 (e.g. 18 s to 22 s) to more clearly visually distinguish the dashed lines from the solid lines.  trajectories: in the top-right of Fig. 10(b), we run inference for our trained model on T2. While the proprioception estimate tracks the general shape of the trajectory well, the optimization, particularly for D y , gets trapped in local minima sometimes leading to periods of higher error. Next, we consider a model trained and evaluated on separated sets of the T5 trajectory. The configuration estimate tracks the ground truth very well as can be seen in the bottom right. Finally, we present a sequence of stills based on camera views and renderings of the soft segment for inference on T3 in Fig. 11.

Conclusion
This work proposed a sensing strategy for soft-bodied robots that relies on multiple magnetic sensors embedded directly in the robot. Thanks to a novel kinematics-aware neural architecture, we can simultaneously use information coming from all the sensors to reliably reconstruct the full robot shape. The decoupling of the kinematics from the learned sensor measurement predictor allows modifications to the placement of the sensors without requiring a re-training of the neural network. The proposed method is agnostic to the used kinematic state parametrization, which we verified in simulations using either PCC or affine curvature models. Extensive experiments with a soft segment showed that a model can be trained on one trajectory type and then be used for inference on a variety of other trajectories in the same workspace. In future work, we will validate the proposed proprioception methodology to execute closed-loop control. We also invite future research studying the optimal placement of sensors in continuum soft robots, where the optimization procedure might take advantage of the gradients provided by our algorithm.

Conflicts of interest
There are no conflicts to declare.
angle y i (q,v) is found by integrating the curvature The rotation to the frame {S v } can then be easily determined with R i iÀ1 ðq; vÞ ¼ R f i ðq; vÞ R y ðq; vÞR T f i ðq; vÞ. After substituting S Á = sin(Á), C Á = cos(Á) for conciseness, we state the homogeneous transformation as We choose to integrate the translational terms in (18) numerically with 101 sample points using the Torchquad 37 implementation of the Simpson's rule, which makes the forward kinematics fully and automatically differentiable.