## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Shape-programmable magnetic soft matter

Edited by John A. Rogers, University of Illinois, Urbana, IL, and approved August 8, 2016 (received for review May 22, 2016)

## Significance

At small scales, shape-programmable magnetic materials have significant potential to achieve mechanical functionalities that are unattainable by traditional miniature machines. Unfortunately, these materials have only been programmed for a small number of specific applications, as previous work can only rely on human intuition to approximate the required magnetization profile and actuating magnetic fields for such materials. Here, we propose a universal programming methodology that can automatically generate the desired magnetization profile and actuating fields for soft materials to achieve new time-varying shapes. The proposed method can enable other researchers to fully capitalize the potential of shape-programming technologies, allowing them to create a wide range of novel soft active surfaces and devices that are critical in robotics, material science, and medicine.

## Abstract

Shape-programmable matter is a class of active materials whose geometry can be controlled to potentially achieve mechanical functionalities beyond those of traditional machines. Among these materials, magnetically actuated matter is particularly promising for achieving complex time-varying shapes at small scale (overall dimensions smaller than 1 cm). However, previous work can only program these materials for limited applications, as they rely solely on human intuition to approximate the required magnetization profile and actuating magnetic fields for their materials. Here, we propose a universal programming methodology that can automatically generate the required magnetization profile and actuating fields for soft matter to achieve new time-varying shapes. The universality of the proposed method can therefore inspire a vast number of miniature soft devices that are critical in robotics, smart engineering surfaces and materials, and biomedical devices. Our proposed method includes theoretical formulations, computational strategies, and fabrication procedures for programming magnetic soft matter. The presented theory and computational method are universal for programming 2D or 3D time-varying shapes, whereas the fabrication technique is generic only for creating planar beams. Based on the proposed programming method, we created a jellyfish-like robot, a spermatozoid-like undulating swimmer, and an artificial cilium that could mimic the complex beating patterns of its biological counterpart.

Shape-programmable matter refers to active materials that can be controlled by heat (1⇓⇓⇓–5), light (6, 7), chemicals (8⇓⇓⇓⇓–13), pressure (14, 15), electric fields (16, 17), or magnetic fields (18⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–33) to generate desired folding or bending. As these materials can reshape their geometries to achieve desired time-varying shapes, they have the potential to create mechanical functionalities beyond those of traditional machines (1, 15). The functionalities of shape-programmable materials are especially appealing for miniature devices whose overall dimensions are smaller than 1 cm as these materials could significantly augment their locomotion and manipulation capabilities. The development of highly functional miniature devices is enticing because, despite having only simple rigid-body motions (34⇓–36) and gripping capabilities (37), existing miniature devices have already been used across a wide range of applications pertaining to microfluidics (38, 39), microfactories (40, 41), bioengineering (42, 43), and health care (35, 44).

Among shape-programmable matter, the magnetically actuated materials are particularly promising for creating complex time-varying shapes at small scales because their control inputs, in the form of magnetic fields, can be specified not only in magnitude but also in their direction and spatial gradients. Furthermore, as they can be fabricated with a continuum magnetization profile,

Although shape-programmable magnetic soft materials have great potential, previous work can only rely on human intuition to approximate the required **m** and actuating fields for these materials to achieve their desired time-varying shapes. As a result, such heuristic and unsystematic methods can only program these materials for a limited number of applications, demonstrating only either simple deformations (21, 24⇓⇓⇓–28) or very specific functionalities (18⇓–20, 22, 23, 29⇓⇓⇓–33).

Here, we present a universal programming methodology that can automatically generate the required *A*). This universal method therefore has the potential to inspire a wide variety of miniature devices that could transform robotics, material science, and biomedicine. The proposed method consists of theoretical formulations, computational strategies, and fabrication procedures. Although the theory and computational method are universal for programming 2D or 3D time-varying shapes (*S2. Programming Materials with 3D Time-Varying Shapes*), our fabrication technique is only generic for making planar beams. Despite the limitations of our current fabrication technique, it is still significantly better than existing techniques, which at most can only create direction-varying

## Programming Methodology

We demonstrate our programming methodology with a large deflecting beam subjected to quasistatic conditions (Fig. 1*B*). For practical considerations, we have also considered

To simplify our discussion, we constrain the cross-sectional area of this beam to be uniform and allow it to bend in a plane. Furthermore, although we have provided the generic discussions for the beam’s boundary conditions in *S1. Boundary Conditions* (see also Fig. S2), here we simplify the boundary conditions to fixed-free (Fig. 1*Bi*). Without any loss in generality, the bending axis of the beam is described by the *z* axis of the global frame shown in Fig. 1*B*.

### Theoretical Formulation.

Following the steps in Fig. 1*A*, we first define the desired deformations along the beam’s length, *s*. Because such deformations can vary with time, *t*, the kinematics can be mathematically represented with the rotational deflections along the beam, *Bi*). After the kinematics are specified, we establish the torque balance equation for an arbitrary infinitesimal element (Fig. 1 *B, ii*), d*s*, at any time, *t*, to be as follows:

The variables *A* represent the applied magnetic torque (per volume), the beam’s bending moment, and its cross-sectional area, respectively. The other variables, *x*- and *y*-axis internal forces within the beam, respectively. Similarly, the force balance equations of the infinitesimal element can be expressed as follows:*x* and *y* axes, respectively. Thus, by using the Euler–Bernoulli equation and substituting Eq. **2** into Eq. **1**, the desired deflections (i.e., required first derivative of bending moment) can be expressed explicitly by the actuating magnetic forces and torques as follows:

The variables *E*, *I*, and *L* represent the Young’s modulus, the second moment of area, and the length of the beam, respectively. Eq. **3** implies that the material’s desired time-varying shapes can be achieved if the magnetic torques and forces can be programmed to balance the desired first derivative of bending moment, across the entire length of the beam at all times. To determine the necessary

The magnetic torque is a function of

### Computational Method.

In contrast to previous magnetic programming studies, we do not use human intuition to speculate the necessary *A*):

The significant benefit of such representation is that Fourier series is inclusive of all possible discrete or continuous mathematical functions, enabling our proposed method to be universal. The vector *T* represents the total time to complete the shape trajectory. The other variables *i* and *j* are integers that range from 0 to *n* and 0 to *m*, respectively. Thus, by substituting Eq. **6** and Eq. **4** into Eq. **3**, we can obtain the following equation:*s* and *t*, with each 2D Fourier coefficients created from the 1D Fourier coefficients in *S6. Additional Discussion*, we will show the mathematical representation of

By following step 3 in Fig. 1*A*, a computational optimization method is then used to determine the optimal values of the 1D Fourier coefficients to satisfy Eq. **7**. This optimization method is implemented by first discretizing the motion of the beam into *p* time frames, that is, *q* segments in length, that is, *q* new equations for each time frame by substituting different values of *s* along the beam into Eq. **7**. By assembling all of the equations across all time frames, there are a total of

Subsequently, the optimal 1D Fourier coefficients in **10** minimizes the difference between the magnetic actuation and the desired first derivative of the bending moment while subjected to the physical constraints of our systems. After this optimization process has been solved numerically by solvers such as genetic algorithm (45) and gradient-based solvers (46), the optimal

### Fabrication Technique.

Based on the magnitude profile of *x*- and *y*-axis components of *A*; see *Materials and Methods* for more details).

## Results

For the first experimental demonstration of our shape-programming methodology, a millimeter-scale beam was programmed to create a shape that resembled a cosine function when it was actuated by a constant magnetic field (Fig. 1 *C i* and *ii*). Next, we programmed a beam to produce a simple sequence of time-varying shapes with 100 discrete time frames. In each time frame, a uniform curvature was held over the beam, gradually increasing between each frame, until the beam curled into a semicircle (Fig. 2*A*). Despite the large number of time-varying shapes, we obtained a simple **7** (Fig. 2 *B* and *C*). The beam was then fabricated and experimentally manipulated to achieve its desired shapes (Fig. 2*D*). Because the required magnetization profile and actuating fields were relatively simple, we extended this concept to simultaneously control multiple beams that have similar motions. By properly configuring several such beams, we were able to reversibly bend them into a “CMU” logo shape (Fig. 2*E* and Movie S1). We further extended this concept by using two similar beams to form the tentacles of a jellyfish-like robot. These tentacles could generate a fast power stroke and a slow recovery stroke for the robot to swim against the slope of an oil–water interface (average speed, 1.8 mm/s; see also Fig. 2*F*, Movie S2, and *S6. Additional Discussion*, for controlling the stroke speeds). The jellyfish-like robot was also steerable, and these steering strategies are discussed in *S3. Steering Strategies* and Fig. S3.

We had also programmed a spermatozoid-like undulating swimmer. To make this swimming gait more biomimetic than previous spermatozoid-like robots (22), we specified the gait to be a propagating traveling wave, with an amplitude that increases linearly from the fixed end to the free end (Fig. 3*A*). Despite the complexity of the gait, our programming method could obtain the necessary *B* and *C*). After fabricating this swimmer, we experimentally show that it could use this gait to swim on an air–water interface (average speed, 11 mm/s; Fig. 3*D* and Movie S3).

Finally, we created an artificial soft cilium that was able to approximate the complex beating pattern of a biological cilium (47). This beating pattern was divided into two strokes—the power and the recovery strokes (Fig. 4*A*). Due to the complexity of this motion, the optimization problem for obtaining **7** to program them. The obtained results for the artificial cilium are shown in Fig. 4 *B–D* and Movie S4, and the key time-varying shapes that we used to closely mimic the complex beating pattern of a biological cilium were shown in Fig. 4*A*. Although other researchers have had some success in creating time-asymmetrical motions for their artificial cilia (50⇓⇓–53), our artificial cilium is the only one on a millimeter scale that can approximate the motions of a biological cilium.

## Discussion

Although the proposed programming methodology is promising, there are several limitations that need to be addressed in future studies. First, our method cannot produce all possible time-varying shapes when **7**, which can be mathematically described as a function of six sets of 2D Fourier series. As all of the 2D Fourier coefficients in this function are created from the lower-dimensional 1D Fourier coefficients of *s* and *t*. This means that we can only program time-varying shapes whose first derivative of bending moment can be expressed in 2D functions, which are representable by the magnetic actuation function. Although a complete analysis to quantify the range of time-varying shapes achievable by our method is beyond the scope of this paper, we have provided a brief discussion on this topic in *S4. Achievable Time-Varying Shapes*. To better understand the range of achievable time-varying shapes, we will formulate a mathematical model to quantify this range in the future. As another future work, we will also increase the range of achievable shapes by developing more powerful electromagnets that can enable the actuating fields to become position variant (*S6. Additional Discussion*).

Second, because several metastable shapes may exist for a given control input, the programmable material may deform into an undesired shape. However, because the selected metastable shape is highly dependent on the previous shape, this limitation can be moderated by using a finer temporal resolution for the shape trajectories. This moderation reduces the deviation between the desired shape and the previous shape, making it easier to guide the material to deform into the desired shape.

Third, due to finite computational power, it is still an open challenge to obtain a global solution from highly nonconvex optimization problems, which have many design variables (46). As a result, we can only rely on numerical techniques such as the two-step optimization approach to moderate this challenge. Although this approach may allow us to obtain more accurate suboptimal solutions, it may also overconstrain the optimization problem unnecessarily and cause a reduction in the original search space (*S5. Discussion for Two-Step Optimization Approach* and Fig. S4). Therefore, we should only implement the two-step optimization approach when we cannot obtain an accurate solution with one optimization process. In view of this challenge, in the future we will also investigate new numerical techniques, allowing us to obtain more accurate solutions for a highly nonconvex optimization problem.

Fourth, our fabrication technique, the two-step molding process shown in Fig. 5 *A–D*, is only capable of tuning the amount of magnetic particles along the in-plane axis of the beams. Because we cannot tune the amount of magnetic particles along all axes of a structure, we are still unable to create an *F*). We will explore new microfabrication processes that will allow us to create smaller-scale structures with desirable

Fifth, although our proposed computational method is not restricted to creating miniature devices, we have yet to use it for macroscale devices. Nevertheless, the procedures for constructing such macroscale devices are discussed in *S6. Additional Discussion*, and the implementation for these procedures will be explored as a future work.

In summary, we have introduced a universal programming methodology that can enable scientists and engineers to magnetically program desired time-varying shapes for soft materials. The method was validated with a simple showcase, and we demonstrated its versatility by creating a reversible CMU logo, a jellyfish-like robot, a spermatozoid-like undulating swimmer, and an artificial cilium. Compared with other shape-programmable materials that may require minutes to induce a shape change (4, 8), our devices can transform into their desired shapes within seconds. We envision that this methodology may enable researchers to develop a wide range of novel soft programmable active surfaces and devices that can find broad applications in robotics, engineering, and biomedicine.

## Materials and Methods

Here, we provide a detailed discussion for our fabrication technique. The required steps to create the desired magnetization profile for a programmable beam are summarized in Fig. 5. The programmable magnetic soft composite material consists of two components: a passive component and an active component that can be stimulated by magnetic excitation. The active component is created by embedding fine neodymium–iron–boron (NdFeB) particles that have an average size of 5 μm (MQFP; Magnequench) into a soft silicone rubber (Ecoflex 00-10; Smooth-on, Inc.). The volume ratio for the NdFeB particles and Ecoflex 00-10 is 0.15:1. The passive component is created by embedding aluminum (Al) powder with an average particle size of 5 μm into the same type of silicone rubber with the same volume ratio. The volume ratio is selected to ensure that the elastic modulus of the active and passive components are identical, allowing the composite to have a uniform elastic modulus. The relationship between the passive component’s volume ratio and its resultant elastic modulus was experimentally characterized (*S7. Matching the Elastic Modulus Properties* and Fig. S5).

To create a nonuniform *A*). The passive component (in liquid form) is poured into the negative mold and is allowed to cure (Fig. 5*B*). Once the passive component is fully cured, a laser cutter is used to cut out a band with nonuniform width (Fig. 5*C*). The active component (in liquid form) is poured into the mold to replace the removed band (Fig. 5*D*). The two components form a composite of a uniform thickness once the active component is cured. Due to the nonuniform width of the band, the distribution of the active components can be patterned. This allows the beam to have an *E* and *Fi*). Thus, by magnetizing the beam when it is sandwiched between the jigs, the desired orientation profile can be obtained after the applied magnetizing field and the jigs are removed (Fig. 5 *F*, *ii*). The NdFeB particles that are embedded within the active components will be saturated by the large uniform magnetizing field (

## S1. Boundary Conditions

To achieve the desired shapes, the boundary conditions of the beams must always be satisfied. Generally, there are two types of boundary conditions: the fixed and free ends of the beam (Fig. S2*A*). At the fixed ends, the desired vertical and rotational deflections must always be constrained to zero. Conversely, the bending moment at the free end is also constrained to zero because it is difficult to apply torques at the free ends. However, because not all time-varying shapes have a zero bending moment at their free end, we introduce a method to overcome this limitation.

Our proposed method is to physically add a small extension, *B*). We can determine the necessary magnetization profile for this extension by including it into the Fourier series representation of

## S2. Programming Materials with 3D Time-Varying Shapes

Although we did not program materials with 3D time-varying shapes experimentally, we provide the theory and computational method here. Similar to the formulations for programming beams, we first specify the desired time-varying displacement fields along the *x*, *y*, and *z* axes of the material coordinates with the variables

The tensor,

After we determine the stress distribution within the materials, we perform a quasistatic analysis. Because it is more desirable to perform the quasistatic analysis in Eulerian description, we convert the second-PK stress into the Cauchy stress in the Eulerian description as follows:*i*th direction, respectively. These magnetic forces and torques (per volume) are functions of

To account for the change in the magnetization profile after the material has deformed, we must map the magnetization profile from its initial undeformed state to its current deformed state. Because *V*. When a deformation occurs, the magnetic moment changes its orientation, and the magnitude of its volume will also be changed:**S4** and **S6**, we obtain the equilibrium equations expressed explicitly with the magnetic torques and forces.

To have a universal approach to solve for these magnetic torques and forces, we will use Fourier series to represent the magnetization profile and magnetic fields again. Here, we specify the magnetization profile, *x*, *y*, and *z* axes of the material. Similarly, all of the components of the actuating fields, that is, *t* [see more details for the independent spatial gradients of *S8. Experimental Setup*]. Mathematically, this means that *x*, *y*, and *z* axes, respectively. By substituting Eqs. **S12**, **S11**, **S10**, **S9**, and **S6** into Eq. **S4**, we can express the six equilibrium equations in terms of these Fourier coefficients.

Similar to the approach presented in *Computational Method*, we determine the optimal Fourier coefficients via an optimization process. This can be achieved by dividing the motion into *p* time frames, and by discretizing the material’s spatial coordinates into *q* points. By substituting *q* different spatial coordinates of the material into the six equilibrium equations for every time frame, we can create 6*p* × *q* equations that are linearly dependent on the product of the Fourier coefficients:

The vectors *p* × *q* equations, respectively. Subsequently, the optimal 3D Fourier coefficients in **10** from the main text. Once these coefficients are identified, the necessary

## S3. Steering Strategies

There are several strategies to steer untethered miniature devices magnetically. We introduce the first steering strategy by using the jellyfish-like robot as an example. To implement this strategy, we intentionally constrain the magnetization profile of the beams to be symmetrical around the *y* axis of the robot’s body frame (Fig. S3*A*). This allows the robot’s net magnetization to be always parallel to its body frame’s *y* axis. The net magnetization of the programmable material,

In addition to the first strategy, we present two other strategies that allow untethered programmable materials to steer in a plane while still achieving their desired shape transformations. The second strategy is to constrain certain motions of the programmable material so that it is easier to steer the device. The last strategy is to include a rigid component that can be used to control the device’s orientation.

The second strategy can be implemented by placing the material on a liquid interface in which the programmable material is constrained by the surface tension of the fluid. As an illustration, Fig. S3*B* shows the *x–z* plane of the programmable material’s body frame. Due to the surface tension of the fluid, the *z*-axis component of the net magnetization cannot create rigid-body torques that affect the orientation of the material. Thus, the alignment of the material on the liquid interface is solely dependent on the body frame’s *x*-axis component of the net magnetization. As a result, the robot’s orientation can be controlled by using an applied *x*-axis component of the net magnetization. This strategy has been implemented to steer the undulating swimmer (Movie S3).

For the last steering strategy, we can control the orientation of the device by programming the magnetization profile of a rigid component. Multiple feasible magnetization profiles may exist, and an example is shown in Fig. S3*C*. In this case, although the net magnetization for the rigid component is zero, this component can still provide a rigid-body torque that can steer the material. By following the body frame assignment in Fig. S3*C*, a rigid-body torque around the *z* axis can be induced on the rigid component when the spatial gradient, *x*–*y* plane. Thus, by controlling the magnitude of the spatial gradient,*z*-axis torque that is induced by the programmable material. Although this spatial gradient will also induce a *z*-axis force for the *x*-axis component of

## S4. Achievable Time-Varying Shapes

Although a complete analysis for determining the range of time-varying shapes achievable by our method is beyond the scope of this paper, we provide a brief discussion here. The proposed method cannot produce all possible time-varying shapes for small-scale materials because the materials have a time-invariant **7** in the main text suggests that there may not exist a feasible solution when a sequence of time-varying shapes requires a large change in their desired profile of *x*- and *y*-axis components of

## S5. Discussion for Two-Step Optimization Approach

Although the two-step optimization approach can be used as a technique to obtain a more accurate solution from a highly nonconvex optimization problem, it may also overconstrain the problem unnecessarily and cause a reduction in the original search space. As a result, the two-step approach should only be implemented when we cannot obtain an accurate solution with one optimization process.

To facilitate an in-depth discussion, we will use the optimization process for the artificial cilium as an example. If the original one-step optimization approach has been used, its level set can be illustrated with Fig. S4. In this level set, its design variables for *x* axis, whereas the *y* axis shows the design variables for the required actuating fields in the power strokes. The isolines in the level set show the contours of the fitness value, that is, the total error quantified by Eq. **10** in the main paper. This error, *E*, can be expressed as the sum of

When the two-step approach is implemented, the first optimization process can be represented by the plot in Fig. S4 *B*, *i*. The *x* axis for this plot represents the design variables for *y* axis represents *B*, *ii*). Once the optimal design variables in the second optimization process has been determined, we obtain the necessary actuating fields in the power stroke. Because the search space of the second optimization process is also much smaller than the original one-step optimization process, this process also has higher chances to jump out from a suboptimal solution. Because the two-step optimization approach has higher chances to jump out from a suboptimal solution compared with the original one-step approach, it is possible to use the two-step approach to obtain a more accurate solution.

Despite the benefits of the two-step optimization, its total search space is much smaller than the original one-step approach. This can be explained by the constrained relationship between *E*. This implies that, although the two-step approach has simplified the original optimization problem, it may also overconstrain it. Thus, the obtained solution from the two-step approach may not be the global solution that contains the lowest

## S6. Additional Discussion

Here, we discuss the possibility of extending the proposed approach to simultaneously determine the magnetization profiles for multiple beams and the effects when the actuating fields can be varied locally in space. We will also discuss how we can change the speed for the shape trajectory. Finally, we will give the mathematical expression for all of the 2D Fourier series listed in Eq. **7** in the main paper.

The programming method can determine the magnetization profiles for multiple beams simultaneously by using a corresponding set of Fourier series for each beam’s magnetization profile. For example, if there are *r* numbers of beams, there will be *r* sets of Fourier series. Thus, for each time frame, we can create *s* across each beam into Eq. **7** in the main text. By assembling all of the equations across all time frames, there will be a total of **10** in the main text, the optimal Fourier coefficients can be determined, thus generating the necessary magnetization profiles for all of the beams.

On the other hand, if *l* regions, there will be *l* number of independent *l* sets of them. However, Eq. **7** in the main text will have to be slightly modified as we have to substitute the corresponding **10** in the main text. The beam will be able to achieve a larger range of time-varying shapes when the actuating fields become position variant. Note that it is easier to create position-variant actuating fields when the materials are in macroscale. This is because a macroscale device will only require a small magnitude of magnetic field spatial gradients to achieve such actuating fields. The steps to program such macroscale devices are similar to programming microscale devices that have a position-variant

Here, we also discuss how we can change the speed of the time-varying shapes. There is an upper speed limit for the programmable material to change its shape. This limit is dictated by either the speed of the electromagnetic coil system that generates the actuating fields or the fundamental natural frequency of the material. In our experiments, it is the speed of the electromagnets that limits our bandwidth to be 25 Hz. Based on this limitation, we have constrained the fastest component of the Fourier series representing

Finally, we will provide the mathematical representation of

## S7. Matching the Elastic Modulus Properties

Because of the embedded metal particles, the elastic modulus of the composite materials is different from that of pure Ecoflex. The embedded aluminum and NdFeB powders were selected to have the same mean particles size of 5 μm. As the volume ratio of the embedded NdFeB powder to Ecoflex in the active component was predetermined, the active component’s elastic modulus was fixed. Therefore, the elastic modulus of the passive component, Ecoflex with embedded aluminum powder, was tuned by changing the volume ratio of the particles to Ecoflex. The elastic modulus of both the passive and active components were evaluated with a tensile testing machine (Instron 5943; Instron, Inc.). Each volume ratio was evaluated with three experiments, and a linear model was fitted to represent the relationship between the elastic modulus and the volume ratio. Based on the fitted model, the necessary volume ratio for the passive component’s elastic modulus to match the active components was determined (Fig. S5). The corresponding mass ratio was obtained by the following:

## S8. Experimental Setup

The magnetic field and its spatial gradients were generated by an electromagnetic coil system that has eight coils, as shown in Fig. S1. The coil system can be controlled to generate the desired magnetic field and its spatial gradient in the workspace with a uniformity above 95% across a ^{3} volume. The mapping between each coil’s electrical current and the resulting actuating fields can be approximated in a linear equation:

## S9. Parameters for Each Case

Here, we provide the parameters used for each experimental showcase, that is, the dimensions of the beams and the number of Fourier series coefficients, or *n* and *m*, respectively. The number of design variables for each optimization process will also be given. For the artificial cilium, because we have to optimize

## Acknowledgments

We acknowledge Carmel Majidi for modeling discussions and Burak Ozdoganlar and Emrullah Korkmaz for their help with fabricating the precision molds for our soft devices. Furthermore, we thank Steven Rich, Zeinab Hosseini-Doust, and Lindsey Hines for their useful suggestions on the paper's language. We also thank members of the NanoRobotics Laboratory at Carnegie Mellon University and the Physical Intelligence Department at the Max Planck Institute for Intelligent Systems for their helpful discussions. This work was partially supported by National Science Foundation National Robotics Initiative Program Grant 1317477.

## Footnotes

↵

^{1}G.Z.L., Z.Y., and X.D. contributed equally to this work.- ↵
^{2}To whom correspondence should be addressed. Email: sitti{at}is.mpg.de.

Author contributions: G.Z.L., Z.Y., X.D., H.M., and M.S. designed research; G.Z.L., Z.Y., X.D., H.M., O.E., and W.H. performed research; G.Z.L., Z.Y., and X.D. analyzed data; G.Z.L. and Z.Y. developed the theory; G.Z.L. and X.D. implemented the computational approach; and G.Z.L., Z.Y., X.D., H.M., W.H., and M.S. wrote the paper.

Conflict of interest statement: A European patent application related to this article has been recently filed by the Max Planck Institute for Intelligent Systems.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1608193113/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵.
- Felton S,
- Tolley M,
- Demaine E,
- Rus D,
- Wood R

- ↵.
- Hawkes E, et al.

- ↵.
- Mohr R, et al.

- ↵
- ↵
- ↵
- ↵.
- Mu J, et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Crivaro A,
- Sheridan R,
- Frecker M,
- Simpson TW,
- Von Lockette P

- ↵
- ↵
- ↵.
- Roche J,
- Von Lockette P,
- Lofland S

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Zhang J,
- Diller E

- ↵.
- Diller E,
- Giltinan J,
- Lum GZ,
- Ye Z,
- Sitti M

*Proceedings of Robotics: Science and Systems*(Berkeley, CA). Available at www.roboticsproceedings.org/rss10/p13.html. Accessed July 16, 2014. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Boyd S,
- Vandenberghe L

- ↵.
- Saldarriaga J,
- Berger JD

- ↵.
- Rozvany GI,
- Olhoff N

- ↵
- ↵.
- Shields AR, et al.

- ↵
- ↵.
- Vilfan M, et al.

- ↵
- ↵.
- Kiral E,
- Eringen AC

- ↵

## Citation Manager Formats

## Sign up for Article Alerts

## Jump to section

- Article
- Abstract
- Programming Methodology
- Results
- Discussion
- Materials and Methods
- S1. Boundary Conditions
- S2. Programming Materials with 3D Time-Varying Shapes
- S3. Steering Strategies
- S4. Achievable Time-Varying Shapes
- S5. Discussion for Two-Step Optimization Approach
- S6. Additional Discussion
- S7. Matching the Elastic Modulus Properties
- S8. Experimental Setup
- S9. Parameters for Each Case
- Acknowledgments
- Footnotes
- References

- Figures & SI
- Info & Metrics