():– Draft: Comparison and Experimental ©The Author(s) 2019 Validation of Predictive Models for Soft, arXiv:1902.00054v1 [cs.RO] 31 Jan 2019 Fiber-Reinforced Actuators Audrey Sedal, Alan Wineman, R. Brent Gillespie, C. David Remy∗ Abstract Successful soft robot modeling approaches appearing in recent literature have been based on a variety of distinct theories, including traditional robotic theory, continuum mechanics, and machine learning. Though specific modeling techniques have been developed for and validated against already realized systems, their strengths and weaknesses have not been explicitly compared against each other. In this paper, we show how three distinct model structures —a lumped-parameter model, a continuum mechanical model, and a neural network— compare in capturing the gross trends and specific features of the force generation of soft robotic actuators. In particular, we study models for Fiber Reinforced Elastomeric Enclosures (FREEs), which are a popular choice of soft actuator and that are used in several soft articulated systems, including soft manipulators, exoskeletons, grippers, and locomoting soft robots. We generated benchmark data by testing eight FREE samples that spanned broad design and kinematic spaces and compared the models on their ability to predict the loading-deformation relationships of these samples. This comparison shows the predictive capabilities of each model on individual actuators and each model’s generalizability across the design space. While the neural net achieved the highest peak performance, the first principles-based models generalized best across all actuator design parameters tested. The results highlight the essential roles of mathematical structure and experimental parameter determination in building high-performing, generalizable soft actuator models with varying effort invested in system identification. 1. Introduction Soft robots make use of elastic behavior in their constituent materials and therefore often encounter physical phenomena that are neglected in rigid robotic theory. Since soft robots are made from deformable materials and deform themselves, such behaviors include unexpected relationships between expanding and constraining structural elements, as well as deformationand direction-dependent nonlinear stiffnesses. Accurate soft robot models need new fundamental assumptions that capture multi-dimensional elastic behavior. Further exacerbating the modeling need is the great variety of available soft system designs and control techniques that may require models. Authors publishing under the soft robotics umbrella make use of a broad set of structural schemes at both the actuator level (e.g. cable robots vs. fluid-driven robots) and the system-wide level (e.g. rigid connection between actuators vs. monolithic soft systems), as well as broad choices of materials, models, and functionalities (Rus and Tolley, 2015). Not only is there a modeling need in soft robotics, but there is a problem of model choice. Successful soft robot modeling approaches appearing in recent literature can be broadly placed into two groups: datadriven strategies and first principle-based strategies. Data-driven strategies have been used on physically realized systems ∗ The authors are with the University of Michigan – Ann Arbor, Department of Mechanical Engineering. 2 - () to characterize dynamic behavior, and, in some cases, to develop control policies. Diverse schemes have been proven on a variety of soft systems, including neural networks (Gillespie et al., 2018; Giorelli et al., 2013, 2015; Thuruthel et al., 2018), genetic algorithms (Giorgio-Serchi et al., 2017), and other regression techniques (Bruder et al., 2018b; Elgeneidy et al., 2018; Veale et al., 2018). Within first principles-based models, two subcategories emerge. Della Santina et al. (2018a) and others (Bruder et al., 2017, 2018b; Della Santina et al., 2018b; Satheeshbabu and Krishnan, 2017; Tondu and Lopez, 2000) show that lumped-parameter models can be used in design and control tasks. Continuum mechanical formulations have recently been demonstrated to be useful in soft system design and control (Coevoet et al., 2017; Connolly et al., 2017; Deimel and Brock, 2016; Goury and Duriez, 2018; Moseley et al., 2016; Polygerinos et al., 2015b; Uppalapati et al., 2018). First principles-based models also study specific behaviors of soft actuators such as interactions between components (Tondu and Lopez, 2000; Wang et al., 2017). Many of the models cited above were developed as part of a larger effort to demonstrate capabilities of specific soft systems. Others perform the complementary task of exploring a design space comprised of many possible systems, and determining the optimal design for a given task. Yet, a broad comparison of model structure and features is missing: why, when, and how data-driven models, lumped-parameter models, and continuum mechanical formulations succeed and fail in soft robotics is not well understood. In this paper, we aim to build this understanding by showing how these three model types compare in capturing the gross behavioral trends and specific features of a popular class of soft actuators. First, we developed three distinct soft actuator models —a lumped-parameter model, a continuum mechanical model, and a neural network— that relate the multidimensional loading and deformation of this actuator. These models differ in mathematical structure, reaching from simple linear equations to integral expressions to sums of weighted functions. The models also differ in how many parameters they require to be identified from experiment, and in the physical meaning (if any) of the parameters. Next, we generated benchmark data by testing eight soft pneumatically driven actuator samples with varied design parameters, across 22,880 kinematic configurations and pressures. Finally, we compared the models on the gross behavioral features that they were able to capture, and on their performance at predicting kinetics across the design space. In particular, we investigated a soft, fiber-reinforced, pneumatic actuator known as the Fiber-Reinforced Elastomeric Enclosure (FREE) (Bishop-Moser and Kota, 2015), or Fiber-Reinforced Soft Actuator (FRSA) (Connolly et al., 2017). The FREE consists of a cylindrical, elastomeric tube whose wall is embedded with 1, 2, or 3 families of fibers wound in helices. The helical pitch of these fibers guides the motion that the tube undergoes when internally pressurized. FREEs are used in several existing systems including exoskeletons (Koller et al., 2016; Park et al., 2014; Polygerinos et al., 2015a; Singh et al., 2018), soft manipulators (McMahan et al., 2006; Pritts and Rahn, 2004), grippers and hands (Deimel and Brock, 2016), vibration isolators (Scarborough et al., 2012), bio-inspired slithering systems (Branyan and Menguc, 2018), and parallel groupings that augment force generation (Robertson et al., 2017). FREEs are a strong candidate for comparative study because they have core mechanical features that are shared in other soft robot architectures. The pressurizeable, nonlinear and hyperelastic wall of the FREE is shared by the PneuNet (Shepherd et al., 2011) and other soft inflatable robot architectures, while the constraint provided by the fiber is structurally similar to fiber constraints in soft cable-driven robots like the octopus-inspired robot created by Laschi et al. (2012). Here, we focus on single fiber-family FREEs because their asymmetric fiber arrangement causes coupling between length change and twisting when the FREE is pressurized (Connolly et al., 2017; Sedal et al., 2018). Together, these key behaviors contrast with the behavior of traditional rigid or series elastic actuators: the FREE design exhibits unexpected relationships between expanding and constraining elements, as well as a nonlinear direction-dependent deformation-loading relationship. FREE models hence demand different governing assumptions than a traditional actuator. Therefore, the FREE serves as a useful investigative platform that captures core problems in soft robotics. In Section 2, we present each of the three models and the performance metric on which they are compared. In Section 3, we describe the actuator samples and experiments. In Section 4, we present the experimental results and error analysis of 3 Γ Unloaded configuration Ri Loaded configuration with applied pressure P>Patm Ro L γ M φ [R, Θ , Z ] T φ F ri ro l P [r, θ, z] T Fig. 1. Face, side, and isometric views of a FREE in an unloaded, un-pressurized reference configuration (top) and a loaded, pressurized configuration (bottom). the models. In Section 5, we compare the gross trends and specific features of each model to those of the data, and compare the performance of the models. 2. Modeling Each model was formulated to relate the FREE’s deformation, design parameters, and loading quantities. The unloaded reference configuration of the FREE is a thick-walled, cylindrical tube with length L, inner radius Ri , and outer radius Ro that is wound with a helical fiber family. The angle Γ between a line tangent to the i of the tube defines h helix and the axis the geometry of the fiber family (see Figure 1). The set of design parameters p̄ = Γ L Ri Ro describes the FREE’s unloaded geometry. For simplicity, we limited the degrees of freedom (DoFs) under consideration to the FREE’s axial elongation and to its end-to-end rotation about the longitudinal axis. This is the intended region of motion of this actuator. The radius is left free to expand or contract. We assumed all other DoFs to be physically constrained. Thus, the current, loaded kinematic state of the FREE is defined by its new length l and its end-to-end twist angle ϕ. The vector of generalized coordinates is h iT hence given by ~q = l ϕ . The forces corresponding to the coordinates in ~q are the axial load F and the axial moment h iT M , given by the vector of generalized forces ~τ = F M . An important consequence of this assumption is that the modeled FREE retains its cylindrical shape under any deformation in our choice of kinematic state ~q. This assumption is valid for cylindrical FREEs which have not undergone a buckling instability. It is also a good approximation for a lightly bent actuator whose bending radius greatly exceeds its cross-sectional radius. Then, for a FREE with a given set of design parameters p̄, the goal of our modelling effort is to characterize the relationship between the kinematic state ~q and the generalized force output ~τ , and to express this relationship as a function of the applied fluid pressure P . Without loss of generality, we chose to formulate this relationship in terms of a force prediction problem and sought to develop and compare different types of models that establish the following functional relationship: ~τ = f (~q, P, p̄). (1) The inverse problem of finding an actuator’s kinematic state ~q given its generalized forces ~τ and design parameters p̄ may be solved numerically, in the case of first principles models, or by re-training, in the case of data-driven models. The three models presented in the following sections differ significantly in their mathematical structure and in the number of free 4 - () parameters that must be identified through experimental data. The physical meaning (if any) of these parameters is also discussed. 2.1. Linear Lumped-Parameter Model The primary assumption of our first model is that fibers of the FREE are inextensible and create a perfect kinematic constraint forming a helix that encloses the elastomeric tube. Considering only the unbuckled case in which the FREE retains its cylindrical shape, we use this assumption to derive an analytic expression for the outer radius ro (Fig. 1) as a function of the state ~q and the design parameters p̄ (Bruder et al., 2018b): ro (~q, p̄) = √ B 2 − l2 . |Φ + ϕ| (2) Here, the fiber length B is given by B = L cos1 Γ and the initial wrapping angle Φ by Φ = RLo tan Γ. Neglecting the wall thickness of the FREE (i.e., assuming that the inner radius ri ≈ ro ), we can employ this expression to compute the volume V of the fluid inside the FREE: V (~q) = πlri2 = π lB 2 − l3 . (Φ + ϕ)2 (3) Taking the partial derivative of this expression with respect to the kinematic state ~q yields the fluid Jacobian JV : JV = h 2 2 ∂V B −3l = π (Φ+ϕ) 2 ∂~q i lB 2 −l3 . −2π (Φ+ϕ) 3 (4) The transpose of this Jacobian allows us to compute the generalized forces ~τf luid that are created by the fluid pressure P (Bruder et al., 2018b): ~τf luid = JTV P. (5) These fluid forces are summed with the elastic forces that result from the deformation of the FREE wall in axial compression, axial twist, and radial compression: ~τ = ~τf luid + ~τwall . (6) The second major assumption of this model is that the FREE wall provides a linear elastic response to deformations in ~q. h i Using ∆~q = l − L ϕ for the deformation along the generalized coordinates, we can compute these elastic forces as: ~τwall = −K∆~q = − " ka kc kc kb # ∆~q. (7) The stiffness matrix K is a positive-definite, symmetric matrix that is found by fitting model parameters ka , kb , and kc to experimental data. The diagonal elements in this matrix, ka and kb , approximately correspond to the lumped stiffness of the wall in axial compression and twist, respectively. Because of the fibers in a FREE, it is necessary to also include off-diagonal elements kc in K. As expressed in Eq. (2), motion in l and ϕ induces a change in the radius of the FREE and hence a radial compression of the wall. The elastic response to this compression is transferred by the fibers back into the generalized torques ~τ . This effect causes an elastic coupling between the twist and axial directions, resulting in the off-diagonal terms in K. 5 We can express the complete lumped parameter model for a FREE by combining these expressions: ~τ = JTV P − K∆~q, (8) which is the special form of Eq. 1 for the linear lumped-parameter model. For this particular model, it is worth noting that the radius ro , the volume V , and the fluid Jacobian may have singular terms when ϕ = −Φ = − RLo tan Γ. This deformation corresponds to a configuration in which the actuator has been rotated such that the fibers are parallel to the tube’s central axis. In this case, the radius and internal volume become ill-defined and the fibers can store any amount of tension or compression in axial directions. In this formulation, the model parameters that must be fit experimentally are the three parameters of lumped stiffness ka , kb , and kc . 2.2. Nonlinear Continuum Model In our second model, we predict the generalized forces ~τ in Equation (1) with a continuum-based, non-linear relationship. Though we continue to assume that the FREE is a cylindrical tube, its wall thickness is no longer neglected and the fiber is considered to be extensible. Extensible fibers cannot kinematically define the tube’s geometry; the radius and internal volume of the FREE in this model depend not only on the generalized coordinates ~q, but also on the internal pressure P . Finally, we replace the assumption of a linear elastic response and replace it with a nonlinear response arising from consideration of the wall’s deformation and occupied volume. To establish a well-defined problem, we assume that: 1. The actuator deforms from a thick-walled tube to a thick-walled tube of different dimensions, 2. the volume occupied by the material of the FREE wall stays constant (i.e. the FREE wall is incompressible, as evidenced for rubbers according to Gent (2012)), and 3. the fiber is perfectly embedded into the elastomer and evenly distributed through it, such that the interaction phenomena are homogenized throughout the FREE wall and the fiber remains locally tangent to the elastomeric tube. These assumptions enable a modified implementation of a continuum mechanical framework proposed by Holzapfel et al. (2000). This specific model implementation has been presented before (Sedal et al., 2018) and is summarized here for completeness. Under the assumptions listed above, the FREE transforms from a tube with initial dimensions {L, Ri , Ro } for length, interior radius, and exterior radius, and no end-to-end twist (ϕ = 0) to a tube with new dimensions {l, ri , ro } and end- to-end twist ϕ (Figure 1). To define this deformation, we track an arbitrarily chosen elemental volume in the FREE wall. ~ = [R Θ Z]T in the FREE’s load-free configuration, and coordinates The elemental volume has location coordinates X T ~ The following functions define the new coordinates ~x = [r θ z] in the current, loaded configuration such that ~x = g(X). ~ of any point X ~ in the unloaded configuration: ~x = g(X) r= s R2 − Ri2 + ri2 , λz ϕ θ =Θ+Z , L z = λz Z, where λz := l . L (9) (10) (11) (12) 6 - () σθ z σrr σθθ σzz σzz σrr σθθ Fig. 2. Arbitrary elemental volume of the FREE wall (marked by a star) with relevant stresses shown. The deformation gradient F, which describes the deformation of the particles in the FREE, is defined as:   F= ∂r ∂R r∂θ ∂R ∂z ∂R ∂r R∂Θ r∂θ R∂Θ ∂z R∂Θ ∂r ∂Z r∂θ ∂Z ∂z ∂Z     = R rλz 0 0 r R 0 0 0   . rΦ L  λz (13) In Holzapfel’s framework, composite materials are modeled by considering the “strain energy," i.e. Helmholtz free energy, stored in the deformed material. The previously listed assumptions allow us to write the strain energy Ψtotal of the FREE wall as a superposition of the free energy from the elastomer and the free energy from the fiber: Ψtotal = Ψisotropic + Ψanisotropic . (14) Here, each free energy depends on the deformation gradient F so as to represent particular material behavior: the fiber is a “standard" fiber that is anisotropic in space (Gent, 2012; Holzapfel et al., 2000) and the elastomer is an isotropic neo-Hookean solid (Ogden, 1997). The neo-Hookean solid model introduces a material parameter C1 and an invariant quantity I1 associated with the elastomer’s isotropy (Spencer, 1971), such that: Ψisotropic = C1 × (I1 (FT F) − 3), 2 (15) The standard fiber model gives the energy associated with the fiber anisotropy by introducing a material parameter C2 and an invariant quantity I4 associated with the fiber stretch (Spencer, 1971): Ψanisotropic = C2 × (I4 (FT F, Γ) − 1)2 . 2 (16) The stresses (i.e. internal forces per area) inside the FREE wall are found by taking the derivative of the total Helmholtz free energy Ψtotal with respect to the deformation gradient F and including a Lagrange multiplier b that accounts for the FREE wall incompressiblity. The expression for three-dimensional stress stored in the deformed FREE wall is then:  σrr  σ =  σθr σzr σrθ σθθ σzθ  σrz  ∂Ψtotal T F − bI. σθz  = ∂F σzz (17) Using the continuum approach enables a system of partial differential equations for the stresses which express the equilibrium of the material under ~τ . In cylindrical coordinates and under the symmetry of the FREE, it is possible to 7 manipulate these equations to develop an expression for b. See the appendix of Sedal et al. (2018) for a detailed derivation. Since the FREE’s radii ri and ro are unconstrained, we use the same equations of equilibrium in Sedal et al. (2018), subject to the boundary condition of the FREE’s internal pressure P , to solve for ri : −P = Z ro ri 1 (σrr − σθθ )dr r (18) Ro2 − Ri2 . λz (19) where ro is found using the material’s incompressibility: ro = s ri2 + Having established ri and ro , we compute the axial force and moment on the FREE by integrating the stresses over the radius of the FREE: F = −2π Z ro σzz rdr + πri2 P (20) ri M = −2π Z ro σθz r2 dr. (21) ri Here, the stresses are expressed in terms of p̄, ~q, and P through the use of Eqs. (17), (14), and (13). This establishes a nonlinear expression: ~τ = f c (~q, P, p̄). (22) In this formulation, the model parameters that must be fit experimentally are the two material parameters C1 and C2 which represent the stiffness of the elastomer and fiber. 2.3. Neural Network Model In our third model, we predict the generalized forces in Eq. (1) with a data-driven approach using a neural network. Our network implementation is based on the previous success of neural networks in modeling the kinetics (Giorelli et al., 2013) and statics (Giorelli et al., 2015) of structurally similar robots. In particular, we implemented the neural network presented by Giorelli et al. (2015) using the inputs and outputs of the FREE statics problem as defined in Eq. (1). We refer to this publication for a more detailed description of the neural network, but summarize our implementation below for completeness. The set of inputs to our neural network consisted of the kinematic state ~q, the internal pressure P , and the design parameters p̄ of each sample. The outputs of the neural network were the components F and M of the generalized forces ~τ . This neural network was implemented as a shallow network with a single fully connected hidden layer containing 6 neurons using a hyperbolic tangent activation function. A schematic of the neural net with inputs xi , weights and biases w, u, o and b, and outputs F, M is shown in Fig. 3. The weights wi,j and uj,k , and biases oj and b1,2 were fit experimentally through training with a back propagation algorithm. With our chosen network topology, 7 inputs, 2 outputs, and 6 neurons on a fully-connected hidden layer, the network had a total of 42 weights wi,j , 12 weights uj,k , 6 biases oj and 2 biases b1,2 . This gives a total of 62 parameters that must be determined experimentally. However, unlike in the other models, these parameters have no physical interpretation. 3. Hardware Experiments To determine the predictive ability of each model, we performed a suite of experiments on a set of eight FREE samples spanning the design space under various loading conditions and imposed kinematic states. 8 - () Inputs xi: q p P Σ w11 x1 wi1 . . . xi u11 o1 . . . w1j wij u1k Σ uj1 Σ Σ ujk oj F b1 M b2 Fig. 3. Schematic of the neural network, with inputs vecq, p̄ and P (xi for i ∈ {1, ..., 7}), hidden layer neurons oj with j ∈ {1, ..., 6} and outputs F and M . Sample Γ (◦ ) L (mm) Ri (mm) Ro (mm) 1 15 90.48 4.77 6.13 2 25 120.52 4.77 6.62 3 36 98.42 4.77 6.74 4 40 90.48 4.77 6.13 5 50 120.40 4.77 6.41 6 62 99.00 4.77 6.36 7 73 128.9 4.77 6.40 8 76 103.22 4.77 6.18 Table 1: Design parameters p̄: fiber orientation Γ (◦ ) and initial dimensions of length L, inner diameter (ID), and wall thickness t of each FREE sample. 3.1. Samples All sample FREEs were made from cotton thread (“Aunt Lydia’s”, Size 10) adhered to rubber tubing. Natural rubber tubing with specified 9.5 ± 0.25 mm inner diameter, 1.6 ±2 × 10− 1 mm thickness (Kent Elastomer) was coated with a thin layer of rubber cement (Elmer’s), resulting in sample thickness between 1.3 and 2 mm. During winding, the fiber angle was prescribed by a 3D-printed template inserted into the rubber tube. Fiber spacing was 1.67 ± 0.25 mm. After winding, a thin layer of liquid latex (TAP) was applied by hand to further secure the fibers on the tube. Samples were cut between 8 and 12 cm in length. We created eight samples with fiber orientations Γ spanning the design space Γ ∈ (0◦ , 90◦ ) in increments of roughly 10 . Fiber orientation of the finished samples was measured through a photograph in three locations and averaged. Standard ◦ deviation of the measured fiber orientation did not exceed 1◦ . Though Γ is the main design variable, our samples also differed somewhat in length and wall thickness. The dimensions of length and thickness were measured with a micrometer three times on each FREE sample and averaged. Standard deviation of sample length and thickness measurements did not exceed 0.10 mm and 0.13 mm respectively. Table 1 shows the fiber orientation and initial dimensions of each sample. After fabrication, the samples were fit using parallel zip ties to the barbed side of 9.5 mm (3/8 in) single barb to “1/8 in NPT” style pneumatic fittings. 3.2. Testing Platform Each sample was fitted into a custom-built testing platform (Figure 4) designed to elongate, twist and pressurize the samples while measuring loading at the tip and photographing the sample’s outer wall. 3.3. Testing Protocol Prior to testing, test bed error was characterized. The linear actuator and rotational servo were tested to ensure position control capabilities within 0.146 mm and 0.35◦ respectively. Cross-talk between the load cells was measured using known forces and torques between 0 and 10 N and 0 and 100 Nmm. Error due to cross-talk on the force sensor was 3.43E-3 N/Nmm, and the same error on the torque sensor was 6E-1 Nmm/N. Parasitic friction on the platform was measured by cycling known tensile loads between 0 and 7 N, and was found not to exceed sensor resolution. 9 Camera Linear Actuator Servo Motor Air Inlet Torque Sensor Load Cell Sample Camera Low-Friction Bushing Fig. 4. Test bed. The sample FREE, shown with blue fibers, is mounted in the center. On the left, the linear actuator and servo motor fix the length and twist. On the right, the torque sensor and load cell measure the loads while the FREE is inflated via the air inlet through a flexible, lightweight tube. Between the air inlet and the sensors is a cylindrical teflon bushing. Top and side camera take low resolution video and high resolution photos. Each FREE sample was mounted into the test bed using custom NPT thread attachments and teflon tape. Aligned markings were placed on the FREE sample and mounting points to observe potential slippage. Then, each sample was tested for all possible combinations of the following kinematic states ~q and internal pressures P : ∆l(mm) = {−5, −4, ..., −1, 0, 1, ..., 4, 5} (23) ∆ϕ(◦ ) = {−120, −110, ..., −20, −10, −1, 1, 10, 20, ..., 110, 120} (24) Pin (V × 10−1 ) = {0, 1, 2, 3, 4, 5, 6, 7} (25) Iterating through these configurations, we first commanded the rotational servo and linear actuator to create a desired end-to-end rotation and axial stretch/compression. As shown by Eqns. (23)-(25), each sample was tested in 286 distinct kinematic states ~q. Then, gauge pressure was set by a voltage signal that corresponded to pressures between 0 and 72.5 kPa. Each sample was inflated to a control signal for 72.5 kPa in eight steps, and deflated back to atmospheric pressure in two additional steps. This resulted in 2,860 configurations tested per sample. After the command pressure was reached, each of these configurations was held for 20 seconds while data (force, torque, and pressure) were collected at 1 Hz. This data was synchronized and averaged to yield a single measurement triplet of ~τ , ~q, and P . Throughout each test, the top camera took time-stamped video at 15 fps to capture any unforeseen events. After the 20 seconds, the side camera photographed the FREE sample wall and then the command for the next configuration was sent. Each photo was later manually classified for whether the FREE sample had buckled in the given configuration or not. All buckled samples were excluded from further analysis. Over all eight samples, loading data and images were collected for 22,880 configurations. After testing was completed, samples were inspected and it was verified that none had been damaged throughout the test. 3.4. Model Comparison Procedure Parameter fitting and model evaluation were done on the separated data sets. We randomly partitioned the data obtained from the un-buckled configurations of each sample into a training set (80% of the data) and test set (20% of the data). Thus, 10 - () there were eight distinct, randomly chosen and randomly ordered training sets, and eight distinct and similarly random test sets. We then created two additional training-test pairs by aggregating the training and test data for the even-numbered samples (Samples 2, 4, 6, and 8) and for the data from all the samples. 3.4.1. Error Metric As an aggregated error metric across a set of n data, we expressed the model error E as the root max max mean square error (RMSE) of force and moment errors normalized by their maximal measured values Fmeas and Mmeas , respectively. This error metric was used to fit model parameters to training data and to determine model performance against test data: v u n  2  i 2 i i i u 1 X Fmeas Mmeas − Mmodel − Fmodel t + E= max max n i=1 Fmeas Mmeas (26) 3.4.2. Model Parameter Identification Model parameters were fit to training data by minimizing the error metric E (Eq. (26)). The fitting operation was undertaken for each of the three models individually on all ten training data sets. Each individually trained model, with its corresponding parameter values, was then used to predict the relationship between loading, pressure, and kinematic state for each test data set. Predicted values were then compared with the measured values to determine model performance using E. For the continuum model, we additionally calculated E on all test sets using model parameters gathered from tests on isolated materials. The three model parameters ka , kb , and kc of the linear lumped model were fit through a constrained minimization of E over a given training data set. The required positive-definiteness of the stiffness matrix K was implemented as a positivity constraint on the eigenvalues of K. The constrained optimization problem was solved with a gradient based interior-point algorithm (fmincon) in Matlab. The algorithm was initialized at order of magnitude estimates of ka = kb = kc = 1. In a similar fashion, we found the two parameters of the nonlinear continuum model, using the positivity constraint C1,2 > 0. The initial estimates for C1 and C2 were 105 and 106 Pa respectively. These are order of magnitude estimates based on previous physical measurements of rubber and cotton fibers (Gent, 2012; Sedal et al., 2018). The constrained optimization problem was again solved using fmincon. Since the constants C1 and C2 represent actual physical properties of the material used for the wall and the fiber, they were also identified via experiments on small material samples. In addition to fitting C1,2 from composite sample measurements as described above, we also evaluated the continuum model’s performance using individual constituent values of C1,2 of the elastomer and fiber respectively from Sedal et al. (2018). The 62 parameters of the neural network were fit with back propagation on normalized input values. We performed training using an unconstrained gradient based method (Levenberg-Marquardt algorithm), implemented with Matlab’s train function. Unlike the other models, the neural network requires the data to be paritioned into 3 sets. Thus, we used 64% for training, 16% for validation, and 20% for testing. Up to 1000 epochs were permitted. The training of the neural network was based on E 2 as objective function, which has the same global minimum as E. 4. Results 4.1. Buckling Out of the 22,880 tested configurations, 12,575 conditions (55.0%) were unbuckled (Figure ??), while 10,305 (45.0%) showed signs of either axial buckling (Figure ??), torsional buckling (Figure ??), or both (Figure ??). Since the geometric assumptions of Section 2 are broken in buckled FREEs, they were not included in subsequent analysis. The remaining configurations are shown in Figure 6. It is not too surprising that the FREEs buckled under several tested conditions: the 11 Fig. 5. Images from experiment of sample configurations: (a) un-buckled, (b) axially buckled, (c) torsionally bucked, and (d) both axially and torsionally buckled. +120° Admissible (Un-Buckled) Measurements All 9 φ -120° -5mm l-L 8 +5mm (a) Sample 1 (Γ = 14◦ ) 7 (b) Sample 2 (Γ = 25◦ ) (c) Sample 3 (Γ = 36◦ ) (d) Sample 4 (Γ = 40◦ ) 6 5 4 3 2 1 0 (e) Sample 5 (Γ = 50◦ ) (f) Sample 6 (Γ = 62◦ ) (g) Sample 7 (Γ = 73◦ ) (h) Sample 8 (Γ = 76◦ ) Fig. 6. Diagrams showing admissible data (from un-buckled trials) to be used for analysis. Data are organized by kinematic state ~ q for all 8 samples. The horizontal axis is the length change l − L in mm, while the vertical axis gives the twist ϕ in degrees. At each of these kinematic states, the FREE did not necessarily stay buckled or un-buckled across the pressure range imposed. The relative proportion of admissible (i.e. un-buckled) measurements at each kinematic state ~ q is reflected by the color of the circle. broad, standardized set of kinematic states and pressures imposed on all samples irrespective of their designed operating range included many configurations outside of the typical range of use of a given FREE. We were not able to identify a clear pattern for when buckling occurred. While there appears to be some tendency for buckling under axial compression and negative end-to-end rotations, this did not hold for all samples: buckling was also observed in axial stretch and positive end-to-end rotation. Furthermore, some samples buckled primarily at high input pressures, while others buckled at low pressures. This lack of a pattern is likely attributed to the presence of multiple different modes of buckling that depend on a sample’s initial geometry, loading, kinematic state and fiber angle (Han et al., 2013; Liu et al., 2014), as well as defects and imperfections in the samples (Lee et al., 2016). A deeper investigation and classification of the observed buckling modes was not the focus of this project. 4.2. Actuator behavior & Model Features To highlight the general behavior of the models, we first present a selection of data: two data sets recorded at specified kinematic states with corresponding comparisons of the modeled axial force F and moment M as functions of input pressure P to the measurements (Figures 7 and 8). We then present all the recorded behavior of a particular FREE sample (Figure 9). Figure 7 shows the data for Sample 6 at ~q = [4mm 80◦ ]T (i.e. under axial extension and positive twist) at 10 pressure values between 1.54 kPa and 63.9 kPa. The models shown here were all trained on the data from Sample 6. 12 - () Figure 7a shows the axial force F arranged as a function of internal pressure P , and Figure 7b shows the moment M as a function of P . Figure 7c shows the same measurements arranged to show the F -M relationship parameterized by P . Similarly, Figure 8 shows the data for Sample 3 at ~q = [−5mm 10◦ ]T (i.e. axial compression and positive twist) and five pressure values between 37.3kPa and 64.4kPa, with all model parameters trained on Sample 3. In the data set of Figure 8, five configurations, all at pressures lower than 37.3kPa, have been excluded from the analysis because the sample was buckled. The same format as Figures 7c and 8c is used for the full set of measurements of Sample 3, shown in Figure 9. Here, the F -M relationship is shown as a vector (i.e. the vector generalized forces ~τ ). The correspondence of the data sets is shown by the pink insert, also present in Figure 8d. From the full behavior of the sample in Figure 9, we can observe general system trends including the relative influence of kinematic state and internal pressure on the magnitude and direction of F and M . The forces and torques produced at un-buckled configurations are characterized by behaviors due to the wall’s elasticity, and behaviors due to the pressure input. Force and moment offsets reflect elastic behavior of a FREE at an imposed kinematic state ~q, at P ≈ Patm . Loading trends reflect how the loads ~τ change as a function of pressure input P . Both of these depend on the initial design parameters p̄. We did not observe significant hysteresis in this experiment. For example, in the kinematic configuration called out in Figures 7 and 8, loading measurements ~τ at similar pressures overlap. Further, the paths taken by the vectors ~τ in Figures 7, 8 and 9 do not appear to change significantly during the pressure ascent or descent. In general, the forces ~τ produced by an actuator are highly dependent on its kinematic state ~q. Throughout the experimental data set, larger magnitudes of ~τ tend to occur at higher pressures P and larger deformations ~q. This is reflected for Sample 3, where the magnitudes of ~τ in the upper right and left quadrants of Sample 3’s kinematic space (Figure 9) are larger than those in the center. Further, the largest magnitudes of ~τ of Sample 3 occur in the upper right quadrant, where the kinematic states impose tension on the fibers. In the lower left quadrant, low forces occur near kinematic states where the FREE was buckled for all the pressures P tested (Figure 6c). These low forces may indicate the onset of the buckling instability, where the fibers are no longer in tension but the wall has not yet buckled. At each kinematic configuration, the direction of ~τmeas can vary with pressure. Comprehensive measurements for all samples at all un-buckled configurations are found in the data packet available as a supplement to this paper. Model comparisons to these data are also shown in Figures 7 and 8. The shapes of the pressure-force, pressure-moment, and force-moment relations for the linear model (shown in blue) are all straight lines. The x− and y− intercepts of each of these straight lines are set by the experimentally determined parameters of the FREE wall stiffness matrix K, while the fluid Jacobian JV determine their slopes as a function of pressure P . The slight curves in the continuum model lines (shown in orange) reflect its geometric and strain-energy nonlinearities. The curve produced by the neural network (shown in green) does not extrapolate smoothly, but shows a hooked shape outside of the region of available data in Figure 8. In this instance, the neuron weights and activation function biases fit to the region of available data did not necessarily extrapolate in a comparable way to the physical models. Model performance varied widely across samples, kinematic states, and input pressures. Though the curves of Figures 7 and 8 demonstrate how each model’s mathematical structure affects its behavior, they should not be used to draw conclusions about which model is most accurate across the data set. A gross performance comparison of the models organized by experimental parameter training set is given in Section 4.3. 4.3. Model Comparison After partitioning the data into training and test sets, we evaluated model performance for every possible training-test pair. The result is 100 error calculations per model, which are depicted on the heat maps of Figure 10. Since the parameters 13 20 0 40 60 Moment (Nmm) Force (N) 1 80 -1 -2 20 0 40 60 80 -20 -40 -3 Pressure (kPa gauge) (b) Pressure (kPa gauge) (a) Force (N) -3 -2 0 -1 1 2 1.54 kPa 17.0 kPa 23.2 kPa 28.8 kPa 36.1 kPa -20 43.8 kPa 47.3 kPa 53.1 kPa -40 63.9 kPa Linear Model Continuum Model Neural Network Moment (Nmm) (c) Force (N) 5 0 20 40 60 -5 Moment (Nmm) Fig. 7. Force and torque models (curves) compared to measurement (black squares) for Sample 6 at ~ q = [4mm 80◦ ]. From top left: (a) T force F by pressure, (b) moment M by pressure, and (c) curves of ~τ = [F M ] parameterized by pressure P . In subfigures (a)-(c), solid lines show the predictions of ~τ corresponding to the pressures P ∈ {1.54, 63.9} in which loading measurements were taken. Dashed lines show a larger modeled range for this sample from P = 0 kPa (upper left end) to P = 70.0 kPa (lower right end). 0 -6 -4 40 60 Linear Model -20 Continuum Model -40 Neural Network Pressure (kPa gauge) (b) Pressure (kPa gauge) (a) -8 20 -2 Force (N) 0 2 Sample 3: -10 37.3 kPa 43.6 kPa 46.8 kPa -20 -30 -40 55.6 kPa 65.4 kPa -50 Moment (Nmm) (c) Moment (Nmm) 0 q= -5mm 10º Force (N) 1 2 3 37.3 kPa -20 -40 65.4 kPa (d) Fig. 8. Force and torque models (curves) compared to measurement (black squares) for Sample 3 at ~ q = [−5mm 10◦ ]. From top left: T (a) force F by pressure, (b) moment M by pressure, and (c) curves of ~τ = [F M ] parameterized by pressure P . In subfigures (a)-(c), solid lines show the predictions of ~τ corresponding to the pressures P ∈ {37.3, 64.4} in which loading measurements were taken. Dashed lines show a larger modeled range for this sample from P = 0 kPa (upper right end) to P = 69.4 kPa (lower left ends). 14 - () Force (N) +120º Moment (Nmm) 0 φ 1 2 3 -20 -40 > 65 kPa 55 - 65kPa 45 - 55 kPa 35 - 45 kPa 25 - 35 kPa 10 - 25 kPa <10 kPa -120 º -5 mm +5 mm l-L Linear Model Continuum Model Mean Error: 34.9% Mean Error: 18.1% Neural Network Error (%) Mean Error: 29.6% 5 6 7 8 2,4, All 6,8 200 6.25 2.19 150 2.32 100 3 4 38.8 50 1 2 Test Sample Set Fig. 9. Raw loading data of Sample 3, shown as vectors ~τ for all pressures across the space of deformation states ~ q . The pink square shows one example of how ~τ changes with internal pressure at ~ q = [−5mm, 10◦ ]T . 192 91.0 1 2 3 4 5 6 7 8 2,4, All 6,8 1 2 3 4 5 6 7 8 2,4, All 6,8 1 2 3 4 5 6 7 8 2,4, All 6,8 0 Training Sample Set Fig. 10. Normalized error, as described in Eq. 26, shown as heat maps, for all possible training-test pairs. Within each heat map, each square’s color designates the error value according to the scale shown in the figure, with darker squares indicating higher error. The horizontal axis indicates which training data were used to fit parameters, while the vertical axis indicates the test set used to evaluate performance. Training and test sets were created individually for Samples 1 through 8, and for two composite sets of Samples 2, 4, 6 and 8 and all the samples, respectively. Maximum and minimum normalized errors for each model are shown in percent in the corresponding cells of the heat maps. Mean normalized error of each model is noted above its heat map. 15 of the continuum model may also be fit through separate experiments on the FREE’s individual component materials, 10 additional error figures are shown in Figure 11. The heat map of the linear model (Figure 10, left) is characterized by two low-error rectangular zones: a larger one at the top right and a smaller one at the bottom left. The low-error rectangular zones show regions where the linear model generalizes. When the stiffness parameters in K are fit to any of Samples 2 through 4, the model extrapolates relatively well in the corresponding test sets, shown by the relatively paler regions in the subplot. Similarly, models trained on Samples 5 through 8 extrapolate well across those test sets. Indeed, the lowest error achieved by the linear model is 2.19% for the training and test sets of Sample 8. This is also the lowest error achieved by any of the models for any training-test pair. The separation of this heat map into regions might be due to the Sample 5 reaching the “magic" fiber orientation (Goriely and Tabor, 2013) γ = 54.7◦ that maximizes its internal volume under the linear model’s assumptions. Both theory and observation show that behavior on either side of this “magic" angle will differ as the internal pressure pushing against the sample walls leads to a volume increase. Along the bottom, top, and left edges of the heat map, bands of higher error appear; these occur near the singularity described in Section 2.1. The value of twist that gives a singular configuration of Sample 1 occurs at ϕ = − RLo tan(Γ) = 122.4◦ ; some kinematic states of Sample 1 in this experiment nearly overlap with the singularity. Higher error then occurs in cases where the training set or the test set of from Sample 1. The “All" test set, which includes Sample 1, also has higher error. Here, the fitted parameters of the stiffness matrix K differed by a maximum factor of 686 across all samples (ka trained from Sample 8 vs. ka trained from Sample 1), but by a maximum factor of 4.7 within the range of Samples 5 through 8 (ka of Sample 8 vs. ka of Sample 5) , and a maximum factor of 3.2 between Samples 2 and 4 (ka of Sample 4 vs. ka of Sample 2). The heat map of the continuum model (Figure 10, center) is characterized by its low variation in error: no regions of the map are particularly bright or dark. Instead, regions of relatively high and low error of the continuum model form horizontal “bands" corresponding with the test set used. For example, the error of Samples 7 and 8 is lower than the minimum error of Sample 2, no matter to which measurements the parameters C1,2 were fit. A particularly dark band occurs when the test set used is Sample 4, but the maximum error of this model (38.8% for the model parameters fit from the Sample 2 training set and tested on the Sample 4 test set) is still substantially lower than the maximum error of the other models. This lower range of error is further reflected in a the lower range of parameters fit across the design space. Here, the parameters fit within the continuum model differed across all samples by a maximum factor of 3.6, (C1 ∈ {0.94 × 105 , 3.39 × 105 } for Sample 8 and C2 ∈ {0.996 × 106 , 1.10 × 106 } for Sample 5). The continuum model was also evaluated with the physical parameters C1,2 of individual constituent materials (i.e. elastomer and fiber) in Sedal et al. (2018). These are shown in Figure 11. The minimum error here (6.63%) is similar to the minimum error achieved in the earlier heat map (6.25%), while the maximum error is slightly lower (29.7% vs 38.8%). Material constants used here (C1 = 5×105 Pa and C2 = 1×106 Pa) and the material constants identified in our experiment were within an order of magnitude. The heat map of the neural network (Figure 10, right) is characterized by its bright, low-error diagonal. Performance of the neural network is especially strong when training and test data from the same sample are used, shown by the lower errors on the diagonal of the heat map. Indeed, the neural network produces the lowest error of any model for the test sets of Samples 1 through 7, though the linear model has the lowest error for Sample 8. Further, the relatively lower errors of Samples 6, 7, and 8 persist in the neural network model, shown by the brighter region of the heat map in the upper right hand quadrant. Elsewhere, the neural network has much higher off-diagonal errors. The neural network heat map in Figure 10 provides evidence of over-fitting when aggregate training sets are used. In the training set of a single sample, test conditions (kinematic states ~q and the input pressure P ) are varied but the design parameters p̄ are the same. As a result, weights and biases fit to p̄ are of value zero for all neural networks trained on a single sample. In aggregate training sets, however, the design parameters p̄ vary according to Table 1, forming a relatively sparse space of design parameters compared to the amount of testing conditions. Here, then, the weights and biases fit to p̄ 16 - () Continuum Model Mean Error: 20.6% All 2,4,6,8 Error (%) 200 8 6.63 7 150 6 5 29.7 100 4 3 50 2 1 0 Fig. 11. Normalized error of the continuum model using material properties C1,2 determined on the individual constituent materials in Sedal et al. (2018). Maximum and minimum normalized error are shown on the test set for which they occurred, and mean error is noted above the figure. are likely to be important. Since these weights and biases have no physical meaning, they may not extrapolate or interpolate appropriately to initial fiber orientations, lengths, and radii outside of the trained values. When trained on Samples 2, 4, 6, and 8, the neural network performed relatively poorly at predicting the behavior of Samples 1, 3 and 5, with 192%, 89.8%, and 145% error respectively. The influence of this effective addition of parameters is also seen by comparing neural network error on Sample 5: trained on Samples 2,4,6, and 8, the neural network gives 145% load prediction error on Sample 5. This is higher than its error on Sample 5 when trained only on Sample 2 (54.0%), only on Sample 4 (56.9%), only on Sample 6 (39.6%), or only on Sample 8 (25.9%). Indeed, error here was lower when we effectively ignored the differences in design parameters than when we tried to fit to them. 5. Discussion Though specific modeling techniques have been validated against realized soft robotic systems in the past, a model comparison across classes of models on a common system has not been made. To address this gap, we developed and compared distinct models that relate the loading and deformation of Fiber-Reinforced Elastomeric Enclosures (FREEs). Three static models were developed and evaluated: a linear lumped-parameter model, a nonlinear continuum mechanical model, and a neural network. We compared predictions of the three models to 12,575 loading measurements that form the broadest (to the authors’ knowledge) data set of FREE loading. The data set is spanning eight varied designs under over hundreds of distinct kinematic deformations and input pressures. Together, these evaluations of model performance enable a comparison of the models for peak performance, generalizability, and system identification effort. We begin by evaluating the peak performance of each model. Intuition may suggest that each model would perform best when its system identification and evaluation use measurements from the same sample — i.e., the diagonal of the heat maps of error in Figure 10. The neural network heat map displays this pattern in a pronounced way, with a brighter diagonal than the linear or continuum models. The average neural network error along the diagonal for the eight samples is 5.51%, while the same average error is 21.4% for the linear model and 16.87% for the continuum model. In most cases along the diagonal, the neural network achieves the lowest error (11.5%, 7.45%, 6.00%, 5.96%, 5.97%2.32% and 2.37% for Samples 1-7 respectively). The linear model achieves the lowest error (2.19%) on Sample 8. Bounded high performance could be especially useful for building models that support control schemes for already-realized soft systems in clearly defined environments. Giorelli et al. (2015) use the neural network for this task, learning inverse statics of their soft tentacle system from measurements taken across all of its possible kinematic configurations. 17 Having shown that peak performance of a model may be confined, we now consider how well the models generalize across the eight samples. A model that generalizes well maintains its performance even as its training and test set vary, resulting in a uniform color across its heat map. The continuum model follows this pattern, having the smallest error range of all models across the training-test pairs (32.6%) , and the lowest average error of 18.1%. In contrast, the linear model has an error range of 88.8% and a mean error of 34.9%, and the neural network has an error range of 189.7% and a mean error of 29.6%. Instead of contrasting diagonals and off-diagonals, the continuum model heat map has horizontal “bands" of similar color, with Samples 6-8 performing best. These results suggest that the continuum model’s nonlinear mathematical structure is crucial to its performance: neither its parameters C1,2 (as noted in Section 4.3) nor its range of error differ as much as those of the other models across training sets. The linear model also generalizes, but does so in rectangular regions, described in greater depth in Section 4.3. This type of linear model has successfully predicted loads for specified parallel combinations of FREEs (Bruder et al., 2018b), but faces limitations due its singularity and the delineation of the training-test pairs for which it performs best. When training and test samples differ, the neural network fails to generalize. This behavior is indicated by the bright diagonal and contrasting darker off-diagonals on the neural network heat map. The neural network’s higher amount of experimentally determined parameters (62) and their lack of physical meaning are the likely cause of this behavior. Since the network is agnostic to any physical assumptions and the parameters carry no particular physical meaning, it may enhance its performance by learning systematic error associated with each individual sample (e.g. fiber irregularity), or capturing phenomena neglected in the simple first principles-based models. These phenomena, however, may not occur again or in the same ways on a different sample. In general, the performance of all models presented here is best for Samples 6-8 for any training set used. The key distinction of these samples from Samples 1-5 is in the winding angle Γ of the helical fiber (Table 1), though the initial dimensions also differ slightly. Improved performance may, then, may be due to the relatively more important role of a high-angle fiber in constraining the radial expansion of the FREE, making it act more like a cylindrical piston and better constraining the radial expansion of the FREE. Or, higher fiber angles may be easier to manufacture consistently. This phenomenon should be the subject of future study. Though our results show that the continuum model’s structure is crucial in its performance, it is not clear which aspect of its nonlinear structure is most important. In contrast to the nonlinear strain energy models presented here, Coevoet et al. (2017) show that finite element material models with nonlinear actuator deformations but linear stress-strain relationships have sufficient accuracy to enable control of a PneuNet-like soft robot, a soft cable robot, and a compliant mechanism. To our knowledge, a continuum model structure with linear material and nonlinear deformation assumptions has not been evaluated systematically for fiber-reinforced actuators; studying how structural aspects of continuum models of FREEs affect their behavior should be the subject of future work. Finally, we compare the system identification effort necessary in each of the models. A model with more experimental parameters is likely to require more dense experimental data, and hence a longer and more exacting data collection scheme, to perform well. In contrast, a model with an appropriate mathematical structure may perform well with fewer experimentally determined parameters and less data collection. The linear model and the continuum model have 3 and 2 experimentally determined parameters respectively, contrasting with the neural network which has 62 experimentally determined parameters. Even when the linear and continuum model fail, they do not reach error figures as high as the maximum errors of the neural network. As indicated by Figures 7 and 8, the physical assumptions on which these models are built give them mathematical structure which enables extrapolation across kinematic states and input pressures. Figures 10 and 11 show that they can also extrapolate across designs. The continuum model has the most flexibility in system identification. This is because the parameters C1,2 have physical meaning that relates to the stiffnesses of a FREE’s constituent materials, rather than its design. As shown by the comparable errors of the continuum model in Figures 10 and 11, a roboticist could fit C1,2 for any existing FREE or from un-assembled constituent materials, and obtain comparable 18 - () model performance as they might on the exact system that they plan to use. For the linear model, generalization occurs in bounded regimes (shown by the brighter rectangles of the linear model heat map in Figure 10): roboticists, then, are able to perform system identification for this model from a limited choice of other, assembled FREE designs. The key way to improve the performance of the linear and continuum models is to incorporate additional or refined physical phenomena. These new phenomena could be specific to an experiment or context, including interactions between our sample and test bed or the effect of defects in the FREE wall. It is possible that for some phenomena to be incorporated into a model, new experimentally determined parameters will need to be added. And, as parameters are added, roboticists will presumably have to trade between model performance and expedient system identification. However, the results in this paper suggest that an appropriate mathematical structure can reduce the need for experimental parameters. The neural network is likely to require the strongest system identification effort. This paper has shown that the same neural network architecture, with the same number of neurons, can at once outperform first principles-based models and under-perform by over-fitting, even when trained on a larger quantity of measurements (i.e. in the composite training set of Samples 2, 4, 6, and 8). As noted in Section 4.3, the neural network has a strong performance when trained on data for a single sample with a dense array of kinematic states and input pressures, but over-fits when trained on multiple FREE samples, which more sparsely occupy the available design space. Roboticists may need to carefully collect training data that are sufficiently dense and that span the entire set of possible inputs, kinematic states, and designs. Without the availability of dense enough data, roboticists may tune a neural network in several ways; a detailed exploration of this is an avenue for future work. Limitations of our study are discussed below. We chose three models —a linear lumped-parameter model, a continuum model, and a neural network— to span a wide space of first principles-based (for the linear and continuum models) and data-driven (for the neural network) techniques. The linear and continuum models are the simplest of their respective classes, and any improved versions of these will inherit their core structure from the models presented here. Conclusions we draw about these classes of models are likely to hold for refined model versions even as performance may improve. In contrast, the neural network is not equally representative of the broad class of data-driven models. It is a popular choice in recent soft robotics modeling work, but its core features are not necessarily inherited by other data-driven techniques. Other data-driven methods relying on fundamentally different mathematical structures may be better at capturing the FREE’s underlying physics. Sünderhauf et al. (2018) note that deep learning techniques have “been shown to learn predictive physics-based models.” Bruder et al. (2018a) show the potential of the Koopman operator in modeling parallel structures of FREEs and Satheeshbabu (2018) uses deep reinforcement learning to model a manipulator made from FREEs. Further, we chose to study FREEs specifically because they share core structural features with both PneuNet-type soft robots and soft cable robots. Despite these similarities, the modeling styles presented here may have different behaviors when applied to other actuator styles. A similar model comparison on non-FREE actuators is an avenue for future work. The static nature of the experimental validation is another limitation. Though viscoelastic behavior was not evident in this experiment, it may still play a role at smaller or larger time-scales than the one presented here. Further, the addition of significant inertial forces to an experiment would clearly change the loading outcomes. A new experiment would be required for a study of the dynamic behavior of FREEs, but all of the models presented here could be adapted to include dynamics. The linear and continuum models here may be extended to include dynamics (Scarborough et al., 2012), and the neural network could be fit on experimental data that includes dynamic phenomena. Having demonstrated the strengths and weakness of these models on predicting the static behavior of FREEs, a dynamic investigation should be the next step. The work presented here provides a broad experimental benchmark enabling comparison of distinct modeling styles found across soft robotics literature. While the modeling approaches presented here have been previously proven on a variety of individual, fully functional soft robotic systems, a comparison on an expansive common data set had not yet been realized. The results shown here confirmed some behaviors already suggested by intuition – e.g. that models with a higher number of experimentally determined parameters may also have the highest peak performance – but also uncovered 19 trends of model performance and physical system behavior across the fiber-reinforced soft actuator design space that are not observable in isolated cases. This broad understanding can help roboticists in various phases of soft system development and validation: our demonstration of the strengths, weaknesses, and failure points of each of these models can guide model choice articulated systems of FREEs or soft actuators similar to FREEs, guide design as roboticists find which schemes best suit their modeling capabilities, and inform the development of a next generation of improved soft actuator models. Acknowledgements This work was funded by the Toyota Research Institute (TRI). We thank our liaison at TRI, Alejandro Castro, for helpful feedback throughout the data gathering and wiritng process. We would also like to thank Prof. Keith Buffinton (Bucknell University), Daniel Bruder, Nils Smit-Anseeuw (University of Michigan) and the members of FlexLab@EPFL for helpful comments, and Daniel Haley for help with experiments. References J. Bishop-Moser and S. Kota. Design and modeling of generalized fiber-reinforced pneumatic soft actuators. IEEE Transactions on Robotics, 31(3):536–545, 2015. C. Branyan and Y. Menguc. Soft snake robots: Investigating the effects of gait parameters on locomotion in complex terrains. In Intelligent Robots and Systems (IROS), 2018 IEEE/RSJ International Conference on, 2018. D. Bruder, A. Sedal, J. Bishop-Moser, S. Kota, and R. Vasudevan. Model based control of fiber reinforced elastofluidic enclosures. In Robotics and Automation (ICRA), 2017 IEEE International Conference on, pages 5539–5544. IEEE, 2017. D. Bruder, C. D. Remy, and R. Vasudevan. Nonlinear system identification of soft robot dynamics using koopman operator theory. arXiv preprint arXiv:1810.06637, 2018a. D. Bruder, A. Sedal, R. Vasudevan, and C. D. Remy. Force generation by parallel combinations of fiber-reinforced fluid-driven actuators. arXiv preprint arXiv:1805.00124, 2018b. E. Coevoet, T. Morales-Bieze, F. Largilliere, Z. Zhang, M. Thieffry, M. Sanz-Lopez, B. Carrez, D. Marchal, O. Goury, J. Dequidt, and C. Duriez. Software toolkit for modeling, simulation, and control of soft robots. Advanced Robotics, 31(22):1208–1224, 2017. F. Connolly, C. J. Walsh, and K. Bertoldi. Automatic design of fiber-reinforced soft actuators for trajectory matching. Proceedings of the National Academy of Sciences, 114(1):51–56, 2017. R. Deimel and O. Brock. A novel type of compliant and underactuated robotic hand for dexterous grasping. The International Journal of Robotics Research, 35(1-3):161–185, 2016. C. Della Santina, R. K. Katzschmann, A. Bicchi, and D. Rus. Dynamic control of soft robots interacting with the environment. 2018a. C. Della Santina, D. Lakatos, A. Bicchi, and A. Albu-Schäffer. Using nonlinear normal modes for execution of efficient cyclic motions in soft robots. arXiv preprint arXiv:1806.08389, 2018b. K. Elgeneidy, N. Lohse, and M. Jackson. Bending angle prediction and control of soft pneumatic actuators with embedded flex sensors–a data-driven approach. Mechatronics, 50:234–247, 2018. A. N. Gent. Engineering with rubber: how to design rubber components. Carl Hanser Verlag GmbH Co KG, 2012. M. T. Gillespie, C. M. Best, E. C. Townsend, D. Wingate, and M. D. Killpack. Learning nonlinear dynamic models of soft robots for model predictive control with neural networks. In 2018 IEEE International Conference on Soft Robotics (RoboSoft), pages 39–45. IEEE, 2018. M. Giorelli, F. Renda, G. Ferri, and C. Laschi. A feed-forward neural network learning the inverse kinetics of a soft cable-driven manipulator moving in three-dimensional space. In Intelligent Robots and Systems (IROS), 2013 IEEE/RSJ International Conference on, pages 5033–5039. IEEE, 2013. M. Giorelli, F. Renda, M. Calisti, A. Arienti, G. Ferri, and C. Laschi. Neural network and jacobian method for solving the inverse statics of a cable-driven soft arm with nonconstant curvature. IEEE Transactions on Robotics, 31(4):823–834, 2015. 20 - () F. Giorgio-Serchi, A. Arienti, F. Corucci, M. Giorelli, and C. Laschi. Hybrid parameter identification of a multi-modal underwater soft robot. Bioinspiration & biomimetics, 12(2):025007, 2017. A. Goriely and M. Tabor. Rotation, inversion and perversion in anisotropic elastic cylindrical tubes and membranes. Proc. R. Soc. A, 469 (2153):20130011, 2013. O. Goury and C. Duriez. Fast, generic, and reliable control and simulation of soft robots using model order reduction. IEEE Transactions on Robotics, (99):1–12, 2018. H.-C. Han, J. K. Chesnutt, J. R. Garcia, Q. Liu, and Q. Wen. Artery buckling: new phenotypes, models, and applications. Annals of biomedical engineering, 41(7):1399–1410, 2013. G. A. Holzapfel, T. C. Gasser, and R. W. Ogden. A new constitutive framework for arterial wall mechanics and a comparative study of material models. Journal of elasticity and the physical science of solids, 61(1-3):1–48, 2000. J. R. Koller, D. H. Gates, D. P. Ferris, and C. D. Remy. ’body-in-the-loop’optimization of assistive robotic devices: A validation study. In Robotics: Science and Systems, 2016. C. Laschi, M. Cianchetti, B. Mazzolai, L. Margheri, M. Follador, and P. Dario. Soft robot arm inspired by the octopus. Advanced Robotics, 26(7):709–727, 2012. A. Lee, F. L. Jiménez, J. Marthelot, J. W. Hutchinson, and P. M. Reis. The geometric role of precisely engineered imperfections on the critical buckling load of spherical elastic shells. Journal of Applied Mechanics, 83(11):111005, 2016. Q. Liu, Q. Wen, M. Mottahedi, and H.-C. Han. Artery buckling analysis using a four-fiber wall model. Journal of biomechanics, 47(11): 2790–2796, 2014. W. McMahan, V. Chitrakaran, M. Csencsits, D. Dawson, I. D. Walker, B. A. Jones, M. Pritts, D. Dienno, M. Grissom, and C. D. Rahn. Field trials and testing of the octarm continuum manipulator. In Robotics and Automation, 2006. ICRA 2006. Proceedings 2006 IEEE International Conference on, pages 2336–2341. IEEE, 2006. P. Moseley, J. M. Florez, H. A. Sonar, G. Agarwal, W. Curtin, and J. Paik. Modeling, design, and development of soft pneumatic actuators with finite element method. Advanced Engineering Materials, 18(6):978–988, 2016. R. W. Ogden. Non-linear elastic deformations. Courier Corporation, 1997. Y.-L. Park, B.-r. Chen, N. O. Pérez-Arancibia, D. Young, L. Stirling, R. J. Wood, E. C. Goldfield, and R. Nagpal. Design and control of a bio-inspired soft wearable robotic device for ankle–foot rehabilitation. Bioinspiration & biomimetics, 9(1):016007, 2014. P. Polygerinos, Z. Wang, K. C. Galloway, R. J. Wood, and C. J. Walsh. Soft robotic glove for combined assistance and at-home rehabilitation. Robotics and Autonomous Systems, 73:135–143, 2015a. P. Polygerinos, Z. Wang, J. T. Overvelde, K. C. Galloway, R. J. Wood, K. Bertoldi, and C. J. Walsh. Modeling of soft fiber-reinforced bending actuators. IEEE Transactions on Robotics, 31(3):778–789, 2015b. M. B. Pritts and C. D. Rahn. Design of an artificial muscle continuum robot. In Robotics and Automation, 2004. Proceedings. ICRA’04. 2004 IEEE International Conference on, volume 5, pages 4742–4746. IEEE, 2004. M. A. Robertson, H. Sadeghi, J. M. Florez, and J. Paik. Soft pneumatic actuator fascicles for high force and reliability. Soft robotics, 4 (1):23–32, 2017. D. Rus and M. T. Tolley. Design, fabrication and control of soft robots. Nature, 521(7553):467–475, 2015. S. Satheeshbabu. Open loop position control of soft continuum arm using deep reinforcement learning. 2018. S. Satheeshbabu and G. Krishnan. Designing systems of fiber reinforced pneumatic actuators using a pseudo-rigid body model. In Intelligent Robots and Systems (IROS), 2017 IEEE/RSJ International Conference on, pages 1201–1206. IEEE, 2017. L. H. Scarborough, C. D. Rahn, and E. C. Smith. Fluidic composite tunable vibration isolators. Journal of Vibration and Acoustics, 134 (1):011010, 2012. A. Sedal, D. Bruder, J. Bishop-Moser, R. Vasudevan, and S. Kota. A continuum model for fiber-reinforced soft robot actuators. Journal of Mechanisms and Robotics, 10(2):024501, 2018. R. F. Shepherd, F. Ilievski, W. Choi, S. A. Morin, A. A. Stokes, A. D. Mazzeo, X. Chen, M. Wang, and G. M. Whitesides. Multigait soft robot. Proceedings of the national academy of sciences, 108(51):20400–20403, 2011. 21 G. Singh, C. Xiao, E. T. Hsiao-Wecksler, and G. Krishnan. Design and analysis of coiled fiber reinforced soft pneumatic actuator. Bioinspiration & biomimetics, 13(3):036010, 2018. A. Spencer. Part iii. theory of invariants. Continuum physics, 1:239–353, 1971. N. Sünderhauf, O. Brock, W. Scheirer, R. Hadsell, D. Fox, J. Leitner, B. Upcroft, P. Abbeel, W. Burgard, M. Milford, and P. Corke. The limits and potentials of deep learning for robotics. The International Journal of Robotics Research, 37(4-5):405–420, 2018. T. G. Thuruthel, E. Falotico, F. Renda, and C. Laschi. Model-based reinforcement learning for closed-loop dynamic control of soft robotic manipulators. IEEE Transactions on Robotics, 2018. B. Tondu and P. Lopez. Modeling and control of mckibben artificial muscle robot actuators. IEEE control systems, 20(2):15–38, 2000. N. K. Uppalapati, G. Singh, and G. Krishnan. Parameter estimation and modeling of a pneumatic continuum manipulator with asymmetric building blocks. In 2018 IEEE International Conference on Soft Robotics (RoboSoft), pages 528–533. IEEE, 2018. A. J. Veale, S. Q. Xie, and I. A. Anderson. Accurate multivariable arbitrary piecewise model regression of mckibben and peano muscle static and damping force behavior. Smart Materials and Structures, 27(10):105048, 2018. Z. Wang, P. Polygerinos, J. T. Overvelde, K. C. Galloway, K. Bertoldi, and C. J. Walsh. Interaction forces of soft fiber reinforced bending actuators. IEEE/ASME Transactions on Mechatronics, 22(2):717–727, 2017.