As previously mentioned, conventional uni-axial mechanical tests on soft materials deform standardized samples over several times their original length (Marechal et al., 2020). However, in more realistic soft robotic applications, especially those driven via magnetic fields, the material may experience only a fraction of these deformations. To improve fidelity between testing for characterization and application, we propose simulating the application-specific actuation method (i.e., magnetically induced forces and torques) within the FEA simulation and solving an optimization problem using an EA to characterize the material. Due to the nature of the EA, the simulation model may be invoked many times, with the potential to significantly increase computational time. Therefore, we developed our FEA model with a balance between accuracy and computational cost to simulate the magneto-mechanical response of our soft robotic designs. This combines analytical mechanical models used to predict nonlinear elastic behavior of polymers with magnetic interaction forces and torques based on a dipole approximation (Petruska and Abbott, 2013; Abbott et al., 2020).

For the mechanical model, we consider soft materials to be homogenous and isotropic. In previous work, it has been shown that compressibility and viscoelasticity may be neglected for lower strains (Steck et al., 2019). In general, to represent the deformed configuration χ of a soft material the following Cauchy tensor is used (Wang et al., 2020): σ = J − 1 ∂ W ∂ F F T , (1) where F  =   ∇ χ is the deformation gradient of the deformed configuration χ , J = d e t  F is the volumetric Jacobian of the deformation and W is the strain energy function that characterizes the specific material. If the material is considered incompressible the Jacobian will be equal to 1. Many different energy functions have been proposed to characterize the behavior of soft hyperelastic materials (Mihai and Goriely, 2017). For this work, we focus on models that have been shown to work well for lower strains (Steck et al., 2019): Neo Hookean, Mooney-Rivlin, Yeoh and Ogden of the third order. These material models, up to the third order, give a wide search space whilst also testing the fitting capability of the framework. The functions considered in the characterization framework are presented in Table 1 in terms of the principal stretches λ 1 , λ 2 , λ 3 . The energy functions each contain a series of constant parameters C_i and α_i that define the specific material. These parameters can be tuned based on experimental observations. These parameters, C_i and α_i, thus constitute the decision space of the optimization algorithm.

The soft material of the FE simulation is then deformed through the interaction with a rigid body, in our case a permanent magnet. The forces f and torques τ on the permanent magnet are induced via a homogeneous magnetic field. These can be modeled by employing the following equations (Jeon et al., 2019; Abbott et al., 2020): f = ∇ ( B ⋅ m ) , (2) τ   =   m   ×   B , (3) where m  is the magnetic moment of the permanent magnet, B  is the magnetic flux density. In a uniform magnetic field, magnetic forces are zero as the field gradient is ∇ B = 0 , leaving only the restoring torque to act on the permanent magnet. In addition, the magnetization of an object is defined as m = v χ ( H ) H where v is the volume of the sample, χ ( H ) is the susceptibility tensor, H  is the applied magnetic field where H = B / μ 0 defines the relationship between the magnetic field and magnetic flux density in which μ 0   is the magnetic permeability of free space. In our case, we can simplify these equations without undermining the validity of the model with the following assumptions. First, since we are using a permanent magnet subject to a small magnetic field, this has no influence on the magnetic moment. Second, since the permanent magnet is physically small, it can be approximated as a point dipole by neglecting its geometry and therefore the susceptibility tensor becomes an identity matrix. As a result, the magnetic moment can be simplified as: m = v B μ 0 = v M , (4) A global coordinate system can be established since the magnetic actuation is the result of an external magnetic field interacting with the magnet. We consider actuation within a 2D plane, therefore, only one component of the magnetic field is non-zero. We define this to be the x -axis according to our reference frame ( B x ≠ 0 and B y , B z = 0 ). The SR was fixed at one end and constrained to bend in only one direction by the applied magnetic field, constraining the magnetic moment to the x y plane. The simplified relationship can be seen in Eq. 5. τ z = − m y B x = − | m | | B | sin θ z , (5) Finally (see Sample Fabrication and Experimental Setup section), body forces due to gravity were integrated in the simulation model, to reconcile with the experimental setup.

The geometry and constraints of the MASR were selected to allow reduction of the model to a 2D plane (Figure 2A), and thus to perform 2D simulations with computational efficiency, by considerably reducing the nodes of the mesh and therefore the complexity of the system. The model so defined can be assumed as plane strain, since the cross-sectional thickness is sufficiently large to consider the depth dimension of the strain tensor equal to zero. For each simulation, geometry of the design was varied using 3 parameters: the length and the width of the robot, and the length of the magnet to match the physical dimensions of the MASR. A quadrilateral mesh was then generated based on the planar geometry (see Figure 2B), and the mesh size was determined after a dependency analysis over the magnet displacement for a fixed magnetic flux density of 15 mT. The magnet was positioned within the robot body with the magnetization direction (north pole) oriented towards the distal end. Therefore, the magnetization direction of the magnet and the externally applied field are perpendicular when the sample is in the resting position, thus maximizing torque (see Figure 2B). The simulation model response may then be evaluated for a range of magnetic fields (Figures 2C, 5B), in line with experimental testing. The model was implemented using COMSOL Multiphysics^®. V5.4 (COMSOL, Sweden) and solved using Multifrontal Massively Parallel Sparse direct Solver. The number of quad mesh elements used for each geometry are summarized in Table 2.

To provide experimental data for the optimization, samples were prepared and tested under a varying magnetic field. Four versions of the design were prepared with the geometrical parameters summarized in Figure 2. One of these geometries (Type 1) was used for training the material characterization through the optimization algorithm, while the others were used for validation (Type 2, 3, 4).

Molds were created for each geometry using a 3D Printer (Ultimaker S5), as shown in Figure 3. Samples were created using two soft elastomers: Dragon Skin™ 10 Medium and Ecoflex™ 00-50 (Smooth-on Inc., United States). For each material, the two-part components were mixed in equal weight using a high vacuum-mixer (Arv-10 from THINKYMIXER, Japan), at a pressure of 20 kPa for 90 s with the centrifuge set at 1,400 rpm to thoroughly mix and remove any air bubbles. A volume of material equal to two times the size of the specimen was slowly injected from the bottom of the mold and residual bubbles and excessive material were expelled using a port at the top of the mold. The material was allowed to cure at room temperature as specified by the manufacturer (5 h for Dragon Skin™ 10 Medium and 3 h for Ecoflex™ 00-50). Once cured, the samples were demolded and the appropriately sized magnets (grade N52, K&J Magnetics Inc., United States) were inserted in the pre-allocated space, with the magnetization direction (north pole) pointing towards the distal end of the MASR.

As shown in Figure 4, samples were tested using a uniform magnetic field, generated using a 1D Helmholtz Coil (DXHC10-200, Dexing Magnet Tech. Co., Ltd., Xiamen, China). During the experiment, the current in the coil was increased in discrete amounts. As a result, the magnetic field was increased proportionally to the current in the coil as described by Eq. 6 (Abbott et al., 2020). B x = ( 4 5 ) ( 3 / 2 ) μ 0 n I R , (6) where n is the number of coil windings, I is the current in the coils and R is the radius of the coils. For each discrete current step, the sample deformation was recorded using a camera (acA 2040-120um, Basler AG, Germany) and post-processed (Figure 4). The resolution of the images was 0.043 px/mm.

To extract image features intrinsically rich in information that could be used to compare the experimental deformation with the simulation model, we developed a post-processing procedure to measure the outer edges of the sample. Bespoke image processing code was developed using MATLAB (Image Processing Toolbox, MathWorks, United States), to allow segmentation and edge extraction where the deflection/deformation is more prominent along the longitudinal length of the sample. The following steps were followed: 1) high contrast images were obtained; allowing the sample edges to be easily extracted; 2) the boundary of the object was extracted, and the corners of the sample were determined from the edges; and 3) after splitting the edge of the sample into different segments, the relevant ones were saved to be used as a target in the characterization process (see Figure 5A).

After the simulated and experimental results were obtained, a metric was needed to compare the two and allow for optimization. We propose the following approach to solve the problem where for each value of the magnetic flux density ( B i ) tested, the predicted simulation model and the recorded experimental data can be compared. First, the distances (or error) d j i between each point j of the simulation contour and corresponding closest point of the experimental sample were evaluated (see Figure 5C). The number of points were in direct correlation with the mesh, such that from geometry Types 1 to 4, the number of points j considered were 80, 110, 70, and 100 respectively. The set of distances created have high dimensionality and thus are impractical to use. Therefore, we used a statistical method to evaluate the fitting errors between the FEA model and the experimental data. The standard error S E i is evaluated for each equilibrium state ( i ) for the magnetic field B i and the chosen model parameters p = [ C 1 … C N , α 1 … α M ] used for the simulation, such as: μ d i ( p , B i ) = 1 N i ∑ j d j i ( p , B i ) , (7) σ d i ( p , B i ) = 1 N i ∑ j ( d j i ( p , B i ) − μ d i ( p , B i ) ) 2 , (8) S E i ( p , B i ) = σ d ( p , B i ) N i , (9) where μ d i and σ d i are, respectively, the mean and the standard deviation of the distances d j i , and N i is the number of points defining the contours of the simulation model for equilibrium state i . After evaluating the standard error for the equilibrium state, we obtained a set of M data (see Figure 5D). Changing parameters for a hyperelastic model (or the model type) may result in a stiffer or softer material compared to the experimental data. This translates into either smaller or larger deflections, making the simulation model drift from the experimental data. To further reduce the dimensionality of this minimization, we decided to evaluate the mean of the standard error ( μ S E ) across all equilibrium test conditions: μ S E ( p ) = 1 M ∑ i S E i ( p , B i ) , (10) The metric in Eq. 10 was used to compare the different material models and as an objective function for the parameter fitting. In addition, the standard deviation ( σ S E ) of the standard error was also used for comparison purposes in the experimental evaluation section (see Experimental Evaluation section): σ S E ( p ) = 1 M ∑ i ( S E i ( p , B i ) − μ S E ( p ) ) 2 , (11)

Here, we defined an optimization problem that allows fitting of the material parameters for a specific analytical model (Table 1). With the goal of minimizing the error by varying the parameters p   of the hyperelastic model in the search space, the following minimization problem was defined: min x μ S E ( p ) subject to  p l ≤ p ≤ p u , (12) where the parameters of the hyperelastic models are subject to upper and lower bounds based on physical constraints. The dimension of p is naturally connected to the order of the model selected. To reduce the search space for p we upper bound the vector element of p to 1 MPa. The objective function involves computing the simulation model using the parameters p and then comparing these with the experimental results using Eq. 10. To solve the defined problem, conventional approaches often rely on the use of simple relations that can be solved using gradient methods such as the Levenberg-Marquardt algorithm (Hartmann et al., 2003; Schumacher et al., 2020). However, due to the complexity of the simulation model, these methods are infeasible for our approach since they fail to navigate over a complex objective landscape. Therefore, we opted to use an EA to perform the optimization.

EAs are, in general, based on some stochastic phenomena. This allows a random walk in the search space to improve the exploration features. Furthermore, these types of algorithms do not require complete information about the problem landscape. These characteristics made this type of algorithm an ideal solution for the problem defined here. In general, EAs are defined as a cyclic process with specific steps (see Figure 6). Initially, a random solution or set of solutions is created, then a two-step process is repeated for several generations (cycles) until the result converges. The first step involves a selection process where the best model parameters are chosen based on the metric (Eq. 10). These solutions are then used to generate new solutions through a stochastic procedure that will replace the solution discarded in the previous step, The stochastic procedure for generating newer solutions is partially based on the information contained in the solution that survived the selection phase (Deb, 2014).

In our framework, solutions that don’t comply with Drucker’s Stability are immediately discarded in the selection phase, specifically, any result that did not converge to a residual error 1 x 10 − 4 within 25 Newton Generation. These results do not represent a valid solution for our fitness function.

Several EA strategies have been proposed, however, one of the most successful to date has been the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) (Hansen and Ostermeier, 2001). This algorithm has shown great potential for solving complex problems compared to other optimization strategies (Caraffini et al., 2013, 2019). This search strategy selects the best solution ν in the population of η individuals for each cycle. It then uses the distribution of the new population to incrementally update a covariance matrix. The principal components of the matrix are subsequently used to direct the exploration in the search space and generate new solutions. By updating the covariance matrix incrementally, the search leads to solutions that show an improvement compared to previous cycles. This population-based strategy is a good compromise between the exploration of the search space and the exploitation of the best solution. For this reason, we decided to use the CMA-ES as a solver for EIMI. Five parameters need to be defined, with the necessary inputs for the algorithm shown in Table 3. The steps described to reach material parameters are shown in the flowchart in Figure 6. The process was implemented in MATLAB (Matlabworks, United States), with a connection to the FE simulation in COMSOL Multiphysics^® V5.4 using the plugin LiveLink™.

The data from the samples chosen for the material characterization are used to define the optimization problem (see Benchmark Functions section) to solve the CMA-ES. The number of evaluations and generations required by the CMA-ES to converge to an optimal solution and obtain a set of model parameters are shown in Table 4. With an increasing number of parameters in the material model, the table displays the growth in the number of generations that were required to converge. This is expected, as they represent a complex objective landscape that is more difficult to explore. The table also shows an outlier, due to the stochastic nature of the EA. The characterization of Ecoflex™ 00-50 using the Ogden model required only 8 generations, despite the search space having high dimensionality with 6 model parameters.

The parameters for the models obtained using EIMI are shown respectively on the top halves of Table 5 for Dragon Skin™ 10 Medium and Table 6 for Ecoflex™ 00-50. The two tables list the mean of the standard error and mean of the standard deviation using Eqs 10, 11 for each type of geometry and model using both our approach and the conventional approach (detailed below).

The results for Dragon Skin™ 10 Medium (see Table 5) show that there is a superior fit achieved for the Type 1 geometry which is the sample used for training. This result can be easily justified from the nature of the process used for characterization. The EA will tend to overfit the parameters for the sample used for training. We can also see that the standard error for the Neo Hookean, Mooney-Rivlin and Yeoh analytical models falls in the 95% confidence interval, making these analytical models similar in performance, as shown in Figure 7A. The Ogden model, however, falls outside of the confidence interval. The number of generations for the Ogden model is higher but could still be insufficient due to a higher number of parameters. To assess the generality of the model characterizations, we evaluated the errors for the other three geometries using the solution obtained by EIMI. For the Type 2 and Type 3 geometries, the standard error for the analytical models falls inside the 95% confidence interval (see Figure 7A). Finally, for the Type 4 geometry, the fitting error of the model parameters is larger. The errors may be caused by different factors, such as in segmentation of the experimental image, or an undesired torsion of the sample as observed during experimental testing (see Sample Type 4 in Figure 8A) (Lloyd et al., 2021) which is not captured within the 2D simulation.

For Ecoflex™ 00-50 (see Table 6), the proposed approach showed similar results to Dragon Skin™ 10 Medium. For the Type 1 geometry, the method showed a tendency of overfitting the models with the number of parameters less than or equal to 3. The modelling errors were also within a confidence interval of 98% (see Figure 7B). Again, the Ogden model has the worst performance for this geometry despite a fast convergence. The optimizer is able to achieve fast convergence for the Ogden model by increasing the dimension of the decision space; however, this results in a more complex landscape that is difficult to explore since more local minima can influence the search. The mean standard error µ S E falls in the 95% confidence interval for all geometries. A slight increase in the values obtained for geometries Type 2 to Type 4 (see Table 6 and Figure 7B). This result may be justified since the increase in the size of the sample points obtained by the segmentation corresponds to an enlargement of the experimental sample dimension.

To further understand the potential of the EIMI approach, we compared it to the conventional approach of fitting material models to experimental uni-axial tensile test data. Here we take advantage of the data and the code provided by Marechal et al. (2020). For each material model analyzed (Table 1), we fit each of the hyperelastic models analyzed using the same tensile data bounded for different engineering strains. More specifically, we fit the models to consider tensile test data for strains of: 20, 50, 100% of the deformation and the full tensile dataset. The obtained parameters and associated errors are shown in the lower part of Tables 5, 6 for Dragon Skin^TM 10 Medium and Ecoflex^TM 00-50 respectively, and in Figure 7. As can be seen in both cases, the conventional approach typically produces larger errors when compared to the EIMI approach. In twenty-one of the twenty-four cases, the overall improvement with EIMI was above 26%.

For the Dragon Skin™ 10 Medium (Table 5), the Yeoh model fit via the conventional method with strains constrained to 50 and 100% of the original length, produced similar results to EIMI; with the latter showing modest improvement of below 2% against the best two models using the conventional approach. However, choosing the most appropriate engineering strain range for fitting purposes is not a trivial matter. The other analytical models (with different strain ranges) using the same dataset failed to get close to the same results, with the EIMI showing at least a 30% improvement. This underlines the difficulty of choosing the most suitable analytical model and the correct stress-strain selection over the tensile data.

Interestingly, for Ecoflex™ 00-50 (Table 6), the best overall performance with the conventional fitting was obtained again using the Yeoh model using the full range (16 times the original length) of strain data for Type 1 and Type 2 geometries and 100% strain for Type 3 and Type 4. Our EIMI showed a performance improvement of 6.8% compared to the best model which in this case was Ogden. For the rest of the cases, EIMI showed an improvement of at least 26%. This result again illustrates how difficult it is to pre-determine the correct range of tensile data to use in order to fit analytical models.

Overall, it is possible to see that the models fit using a conventional approach over different ranges of stress-strain data failed to converge to similar results. The box plots in Figure 7 show a higher variance between conventional models fit using the same tensile range. On the other hand, all the models obtained using the EIMI converge to similar performance values, validating the procedure. Thus, the proposed EIMI method showed that it is capable of dynamically fitting the models without any previous knowledge of the strain involved.

The performance of the different models can be clearly seen in Figure 8 where we can qualitatively observe that our proposed method shows improved performance. In both pictures, the models obtained using EIMI are compared with the models trained using tensile data with a strain of up to 100%. In general, it is possible to see that EIMI is able to produce a better result regardless of the complexity of the model, with similar responses across the different models. Whereas the model fitted using the conventional approach drifts considerably from the experimental data. For Dragon Skin^TM 10 (see Figure 8A), it is possible to see that the Neo Hookean and the Yeoh models underestimated the stiffness, while the Mooney-Rivlin model overestimated the stiffness. For Ecoflex^TM 00-50 (Figure 8B) all 3 models underestimated the stiffness.

To understand how performance varies with actuation, we analyzed EIMI and the best models found using the conventional approach for each of the magnetic field values tested. More specifically we analyzed the SE (see Eq. 9) at each magnetic field value B i recorded for each type of geometry, as shown in Figure 9. This allowed us to understand how the models drift from the experimental sample with increasing magnetic torque.

The performance of the Dragon Skin™ 10 Medium is shown in Figure 9A. For geometry Types 1, 2 and 3, all the models obtained using our EIMI approach produce a smaller error than the best performer (Yeoh) from the conventional approach. Furthermore, the error is constant with the increasing magnetic field, while for the other models it increases with the magnetic field. For Type 4, we see that the conventional Yeoh and Neo Hookean models have a smaller error across all field values. This was due to an observed twisting in the sample that reduced the bending of the beam, which can be interpreted as a stiffer material in the 2D plane (see Supplementary Video S1). This is why we see the errors obtained using the EIMI approach steadily increasing for the Type 4 geometry.

For Ecoflex™ 00-50 (see Figure 9B), we can see that for the sample used for training (Type 1) the models obtained using EIMI outperform all other models. For Types 1 and 2 samples, the Yeoh and Mooney-Rivlin models fitted over 100% strain of tensile data showed slightly worse performance compared to the EIMI and have a similar trend in the SE with the increasing magnetic field. However, these models outperform the EIMI for higher values of the magnetic field in Type 4, again likely due to the twisting of the sample observed experimentally (see Supplementary Video S1).

Overall, our proposed approach is able to maintain a steady error (Figure 9) with increasing magnetic field showing that it is able to dynamically follow the deflection of the robot.