A dynamic model is developed for small-scale robots with multiple high-frequency actuated compliant elastic legs and a rigid body. The motion of the small-scale robots results from dual-direction motion of piezoelectric actuators attached to the legs, with impact dynamics increasing robot locomotion complexity. A dynamic model is developed to describe the small-scale robot motion in the presence of variable properties of the underlying terrain. The dynamic model is derived from beam theory with appropriate boundary and loading conditions and considers each robot leg as a continuous structure moving in two directions. Robot body motion is modeled in up to five degrees-of-freedom (DOF) using a rigid body approximation for the central robot chassis. Individual modes of the resulting multimode robot are treated as second-order linear systems. The dynamic model is tested with two different centimeter-scale robot prototypes having an analogous actuation scheme to millimeter-scale microrobots. In accounting for the interaction between the robot and ground, a dynamic model using the first two modes of each leg shows good agreement with experimental results for the centimeter-scale prototypes, in terms of both magnitude and the trends in robot locomotion with respect to actuation conditions.

## Introduction

This work examines the interaction of structural and contact dynamics in locomotion of small (<10 cm), elastic terrestrial robots operated near resonance. Many approaches have been applied to the general study of robot walking dynamics, but small-scale and/or elastic behaviors can add complexity. For large-scale robots and structures, two approaches to dynamic modeling are extensively used: rigid body models and elastic models. Most dynamic models of large-scale robots use lumped mass, spring, and damper approximations [1,2], though continuous elastic models have also been used for some structures with compliant behavior, typically robot arms described in modal form [3,4].

While large-scale robots are more likely to be actuated with motors at distinct joint locations, small-scale robots often feature compliance throughout or over large portions of their structure. Additionally, locomotion is influenced by light damping compared to larger systems, large fabrication uncertainties relative to the size of structural features, and relatively high resonant frequencies [5,6]. This is especially prevalent in robots based on smart materials and/or micromachined structures. Integrating the resulting vibration behavior with other dynamic interactions encountered by small-robots, such as foot–terrain impact, has yet to be studied in detail.

Understanding the dynamics of walking microrobots is important for design, estimation, and control of efficient locomotion [7]. However, what work has been done on dynamic, small-scale robot locomotion has varied dramatically based on robot design, actuation principles, and operating environment. For example, a previous small-scale (millimeter- to centimeter-scale) robot study showed the importance of dynamic analysis for a type of walking microrobot with multiple pairs of legs, but the study was performed under assumptions of ideal foot–ground contact with excitation from an external vibration field [8]. Other works have measured dynamic performance without detailed analysis [9]. Works focusing on small-scale effects in microrobots using an analytical framework have addressed more limited types of gait or operation, including: (1) stick–slip walking on a smooth surface [10]; (2) multi-axis, multilayer elastic leg dynamic operation in air [11,12]; and (3) ground contact of a single leg [13]. Studies of small robot dynamics at a full-robot level have included finite element analysis [14] and use of lumped parameter models dictating interaction between legs and body [5,6], including with some foot–terrain interaction [15]. A summary of the above combinations of robot and ground contact dynamic modeling techniques is compared in Table 1 for various walking robots with sizes ranging from tens of centimeter to tens of microns.

Unlike these existing works on dynamic locomotion of small-scale walking robots, this paper focuses on the problem of how to effectively account for elastic behavior in miniature robot legs with ground contact. Specifically, we identify an appropriate model for interaction between an elastic robot leg and contact dynamics from ground impact, and highlight some notable resulting phenomena that have influence on the effectiveness or speed of robot locomotion. The model is then verified experimentally using two different centimeter-scale prototypes.

The prototypes used in this study employ piezoelectric actuators bonded to three-dimensional (3D)-printed bodies. Their geometries are designed with legs that have coupled in-plane and out-of-plane dynamics, range of motion, and operating frequencies representative of even smaller, millimeter-scale piezoelectric microrobots made with silicon micromachining [11]. In addition, the two prototypes were fabricated at different size scales, where the larger robot has dimensions of 4 × 7 × 1 cm3, and the smaller robot has dimensions of 1.6 × 3.0 × 0.4 cm3. These small-scale robots have a combination of rigid and compliant structures as well as significant influence of ground contact behavior, resulting in dynamics which demonstrate greater complexity than has been analyzed in prior walking robots of similar size. The proposed model captures motion across multiple legs having multi-axis modal behavior, building upon prior studies of small-scale appendage and single-foot ground contact [13,15]. Some preliminary results have been presented in the authors' previous conference publication [16]. This paper further extends our model with better damping characterization, expanded results for vertical and axial dynamics, and additional results from a smaller prototype to examine scaling effects. The experimental behaviors of the centimeter-scale 3D-printed prototypes are compared with simulations from the dynamic model. Key results include predictions of the direction and amplitude of robot motion under different actuation signal conditions, with more detailed measurements and simulations of vertical motion of the legs. Sensitivity of robot motion over different robot parameters is also analyzed to provide a basis for the design of robot locomotion and parameter optimization.

## Model Description

### General Robot Architecture.

Robots studied in this paper include small-scale (<10 cm) walking robots with multiple pairs of elastic legs connected by a rigid body. Figure 1 shows a schematic view of a robot with 2n legs, each of which may be modeled using stiffness and compliance matrices. The motion of each foot in the y–z plane is examined (in the directions of robot forward and out-of-plane motion from ground). Motion in these 2DOF under proper operating conditions may produce locomotion against friction and gravity forces acting on the robot. In Fig. 1, parameters kj,i and bj,i represent the spring constant and damping coefficient of the ith mode of the jth leg (where i = 1, 2 corresponds to the y- and z-directions, respectively). In addition, each leg is taken to have a known nominal gain with respect to voltage and ratio between motions in y- and z-directions for relevant vibration modes.

This general architecture is inspired by a type of millimeter-scale, piezoelectric and polymer thin-film robot fabricated with micromachining, shown in Fig. 2(a), with a rigid body and multiple pairs of legs on both sides of the body. The legs are made of parylene and silicon microstructures actuated by thin-film lead zirconate titanate (PZT) [11]. Figure 2(b) shows legs designed with 2DOFs, vertical and in-plane, that are connected in series to produce approximately elliptical motion from each appendage. To date, these robots have only been successfully actuated within the silicon chips on which they are built [11]. However, they have shown relatively large ratios of actuation force to the amount of electrical energy required to drive the piezoelectric actuator, making them promising for future autonomous locomotion.

The same robot modal behavior may be realized through other robot designs, which can allow for rapid prototyping and evaluation of dynamics over more varied surfaces and operating ranges. The current study is performed using two different centimeter-scale six-legged prototypes, made with a 3D-printed polylactic acid (PLA) frame bonded with piezoelectric ceramic strips as actuators (Fig. 3). The asymmetric leg actuation of both robots can be modeled with coupled modal dynamics in two directions, the model framework inspired by the millimeter-scale robots [11].

### Leg Model

#### General Description.

Actuation from embedded piezoelectric actuators operating near resonance is intended to generate leg displacements in pseudo-elliptical trajectories, driving the robot forward as contact is made with ground. The current analysis assumes that equal numbers of legs are being actuated in and out of phase at the operating frequency, i.e., exchanging between sets of three legs in a simple hexapod gait.

In the centimeter-scale prototypes developed in this study, the compliant leg takes the form of a misaligned or an asymmetric cantilever beam. Uniform deformation of piezoelectric cantilever microactuators has been studied in some detail [25,26], and beam dynamics are compiled into a conventional modal format in this analysis.

#### Leg Dynamics.

The uniformly distributed piezoelectric actuation force within the beam and point contact force acting on a robot foot, as functions of time, cause the time response of the beam. The analytical solution of the forced response of such a beam, in the vertical direction, is given in a previous study (Eqs. (8.110)–(8.116) in Ref. [27]). The effective force (Qi) acting on the ith mode could be derived as
$Qi=∫0Lf(x)wi(x)dxρA∫0Lwi2(x)dx$
(1)
in which f(x) is the force as a function of position, x, along the cantilever beam; $ρ$ is the density of the beam; and A is the cross-sectional area. With the effective mass for point impulse force ($mi,imp$) and uniformly distributed pulse force ($mi,pls$), the final form of the impulse ($Qi,imp$) and pulse ($Qi,pls$) effective force is
$Qi,imp=wi(L)ρA∫0Lwi2(x)dxfimpδ(t−timp)=1mi,impfimpδ(t−timp)$
(2)
$Qi,pls=∫0Lwi(x)dxρA∫0Lwi2(x)dxFpls=1mi,plsFpls$
(3)

in which $wi$ is the displacement function of ith mode; $ρ$ is the area density of cantilever beam; $A$ is the cross section area of the beam; and $fimp$ and $fpls$ are the impulse force acting on the end of beam and the equivalent piezoelectric force acting on the entire beam.

For derivation of the leg model, modes for each robot leg and foot system are assumed to be independent from each other due to the assumption of a rigid central robot body. Leg models are derived from standard multilayer beam analysis [25] such that a state space form can present the system behavior as follows:
$[X˙j1⋮Xjn]=[Aj100⋱00Ajn][Xj1⋮Xjn]+[Bj1⋮Bjn][fimpδ(t−timp)fpls]$
(4)
$zj=[10…10][Xj1⋮Xjn]$
(5)
(6)

in which j is the leg number; kji and bji are the spring constant and damping coefficient, normalized by effective mass, with respect to each mode for leg j; zj is the jth foot displacement in the vertical direction; $xji$ and $xji˙$ are the displacement and velocity for the ith mode of the jth leg; fimp is the impulse force acting on the tip of leg; fpls is the pulse force acting uniformly on the entire leg; $δ(t−timp)$ is the Dirac delta function representing when impulse force occurs ($timp$); and $mi,imp$ and $mi,pls$ are the effective masses for the ith mode for the impulse and piezoelectric forces. In simulation, the state space form is discretized. Gravitational force on the leg and the internal piezoelectric actuation force are treated as uniformly distributed, while ground contact force and the weight of the robot's foot are treated as acting on a single point at the tip of the leg.

The actuation architecture is confirmed with finite element analysis in comsol multiphysics, as shown in Fig. 4. To perform multi-axis modeling, it is assumed that: (1) the motion in z- and y-directions is coupled, which means their systematic parameters (modal frequency and damping) are identical and (2) the actuation force in z- and y-direction is different, but related by a constant proportionality coefficient.

#### Contact Model: Friction and Restitution.

The contact dynamics of microrobotic appendages at the current size scale was previously studied in Ref. [28]. A continuous contact model is chosen because it is a force-based model for use with the leg dynamics from Eq. (4). To relate contact behavior with measurable velocity data, the robot is modeled with Newton's restitution model, in which the coefficient of restitution is related to the change of velocity before and after contact.

A basic modeling method of restitution and friction is applied. Both forces are assumed to happen during a short period of time. The contact between a foot and ground is modeled simply by the coefficient of restitution ($cr$) from Newton's model and the coefficient of friction ($cf$) from the dry friction model [28]. Considering the velocity in y- and z-direction before ($vyo$, $vzo$) and after contact ($vyf$, $vzf$), behavior during contact is modeled as
$vzf=−vzocr$
(7)
(8)

It is assumed that the ground has infinite mass and negligible compliance, and that the friction does not influence vertical motion. The lateral change in velocity (Eq. (8)) is derived by treating the effective impulse force implied by Eq. (7) as the normal force in friction calculation. The coefficients of certain sample ground surfaces will be identified for simulation in the Model Validation section.

### Body Model

#### Five Degrees-of-Freedom Versus Three Degrees-of-Freedom Interactions.

The robot body could be considered as a rigid body with 6DOFs. The standard moments of inertia in three directions (Ix, Iy, Iz) quantify rotational inertia. The forces acting on the body arise from the spring and damper equivalents in the leg/foot model, acting on the connection points between the legs and the body. Since six legs are symmetrically connected to the body, their forces in two (z and y) directions are considered, but lateral (x-axis) forces are omitted from this analysis.

First, the distances from the body's center of mass to the legs are given in vector form (x, y, z). When calculated for the small-scale robots, this takes the form
(9)
in which lx,i, ly,i, and lz,i are the x-, y-, and z-direction distances from the exterior 2n legs to the body center of mass (the geometry is considered symmetric, which means the two middle legs have y-directional distance equal to zero). The total forces transmitted from the legs for the centimeter-scale robot prototypes are enumerated by vectors
(10)
with Fi,y and Fi,z being the forces from the ith leg acting on the body in the y- and z-directions, respectively, for an n-legged robot. Total moments and forces acting on the body are simply
(11)
and
$Fx=0, Fy=∑Fyl, Fz=∑Fyl−mg−cdvz2$
(12)

Here, Fx, Fy, and Fz are the summed forces acting on the body in x-, y-, and z-directions; m is the mass of the body and g is the gravitational constant. $cdvz2$ is an air viscosity term, which is small, influenced by the air viscosity and the body z-direction velocity. Fx is assumed to be zero since the x-directional translational motion is relatively small from both calculation and experiments. Furthermore, the tilting in the x- and y-direction only has small displacement compared to the overall motion of robots. However, the randomness and uncertainty in x- and y-directional tilting have large influence on the forward motion of the robot since they would change the moments of foot–terrain interaction and when individual leg impacts. Therefore, a 3DOF interaction between body and legs is sufficient to describe the robot at small scale.

#### Relative Motion.

The relative motion between a foot and body, foot and ground, and body and ground dictates the effects of foot impact with ground. The body motion relative to ground is important to reflect body states, which is measured experimentally over a long duration of time. The foot motion relative to ground is necessary for contact performance calculations. The foot motion relative to the body is applied to calculate the force between leg and body. During simulation, for example, once the foot motion relative to the ground is calculated, the foot motion relative to the body is needed to find the influence of foot motion on the body.

Assuming body motion and foot motion relative to ground are known from simulation, the foot motion relative to robot body is determined by the following equations:
$zf=zfb−(zb±xb sin (ϑb,y)±yb sin(θb,x))$
(13)
$vz,f=vz,fb−(vz,b±vx,b cos (ϑb,y)vϑ,by±vyb cos(θbx)vθ,bx)$
(14)
$yf=yfb−(yb±x sin (ϑb,z)±y sin(θb,x))$
(15)
$vy,f=vy,fb−(vy,b±x cos (ϑb,z)v,ϑbz±y cos(θb,x)vθ,bx)$
(16)

in which $xb$, $yb$, $zb$, $vx,b$, $vy,b$, and $vz,b$ are the body displacement and velocity in x-, y-, and z-directions relative to ground; $θbx$, $θby$, $θbz$, $vθ,bx$, $vθ,by$, and $vθ,bz$ are the body rotation and angular velocity in x-, y-, and z-directions relative to ground; $zfb$, $vz,fb$, $yfb$, and $vy,fb$ are the z-directional foot displacement and velocity and y-directional displacement and velocity relative to ground; $zf$, $vz,f$, $yf$, and $vy,f$ are the z-directional and y-directional foot displacement and velocity, relative to body; and x, y, and z are the perpendicular distances between the body mass center and foot location, which is different for each foot. The $±$ signs in the equations are determined by the positive rotational direction of each axis and the location of each foot.

If the body motion relative to ground and the foot motion relative to body were given, the calculation of foot motion relative to the ground would be the backward calculation of the equations shown above. In simulation, the model is implemented using a for loop. Inside each loop, each foot's states determine whether that foot is in contact with ground. If so, the contact simulation is active; otherwise, the simulation mimics the free in-air motion of robot.

## Model Validations

### Parameters: Material Properties, Friction, Restitution; Identified Parameters

#### Detailed Prototype Design.

Both centimeter-scale prototypes are designed with six legs connected to a rigid body; however, each prototype employs a different design for leg geometry. The large centimeter-scale prototype (referred to as the 80 mm prototype in all later sections) is designed with single-beam legs actuated with bimorph piezoelectric ceramic strips, where each beam has a foot at the tip. The piezoelectric ceramic actuators are intentionally misaligned with the PLA leg frame (Fig. 3), which generates a bending moment in the leg with vertical and lateral components. The small centimeter-scale prototype (referred to as the 30 mm prototype in all later sections) uses a different leg design, where two beams are connected with a single foot at the tip. Only one of the beams is bonded with a unimorph piezoelectric ceramic strip, again causing foot motion in two directions. The major characteristics of the 3D-printed centimeter-scale robot prototypes, in comparison to a millimeter-scale version, are described in Table 2. The dynamic model for the millimeter-scale version will be modified from the work presented in this paper. Major parameters of the millimeter-scale version are listed to show the scale difference compared to the prototypes.

Vertical motion of the legs lifts the robot body against gravity, while the in-plane motion is the key feature to move the robot forward. The actuation force tilts the legs about their neutral axis to generate coupled dual-directional motion. The coefficient between the motion amplitudes in the two directions is determined experimentally. For example, the 80 mm prototype with 60 V input has a vertical deflection of about 166 μm near resonance. As shown in Fig. 4, the FEA confirms this actuation behavior for both prototypes. The FEA gives the leg deformation within the same range of measurement. However, due to limited resolution when fabricating prototypes, the FEA is unable to predict the exact dynamics of each individual leg. Therefore, empirically determined parameters are used to model the resonant modes of each leg.

In the full robots, a simple tripod gait is produced by actuating legs 1, 3, and 5 with the same voltage signal while the other three are actuated with the signal out of phase. For the 80 mm prototype, the actuation signal is nominally a 60 V peak-to-peak square wave with 0 V offset, labeled by its 30 V amplitude in later discussion. For the 30 mm prototype, the actuation signal is nominally a 30 V peak-to-peak square wave with no offset.

#### Prototype Fabrication.

The robot prototypes were fabricated using a combination of 3D-printed parts and off-the-shelf components. The frames of the prototypes were 3D-printed with PLA. Next, six piezoelectric ceramic strips (bimorph for the 80 mm and unimorph for the 30 mm prototype) were bonded to the robot frames with epoxy. Additional silver epoxy was used to bond wires to the piezoelectric ceramic strips for actuation. The finished prototypes are shown in Fig. 3.

#### Leg Parameter Identification.

Leg parameters are determined empirically for both 80 mm and 30 mm prototypes. The material properties of the 3D-printed legs are experimentally identified from frequency sweep testing without ground contact. This test is used to identify parameters such as elastic modulus and density of the legs as fabricated. The prototype is fixed with three different boundary conditions as follows: (1) held from both ends on prototype body, (2) held from single end of body, and (3) held at middle of the prototype body. A laser Doppler vibrometer (LDV) and labview are used to record the frequency response of each foot and two points on the body in vertical direction. Experiments are repeated three times.

Figure 5 shows the body vertical motion and two different leg motions under constant boundary condition. Even with some coupling effects between the body and leg, the robot body has relatively little motion in the vertical direction at low frequencies, including through the first three modes of the legs. This allows the assumption of rigid body motion of the prototype body to be used in modeling. The frequency measurement of the same leg under different holding conditions (boundary conditions) appears to have a less than 5% difference in frequency and less than 10% difference in resonance magnitude. With these two holding conditions showing similar frequency responses, this is taken to approximate the free behavior of the robot during locomotion.

The leg parameters of the 80 mm prototype (which are found experimentally) identify the dominant resonant frequencies for six legs (Table 3), one arising in each leg. Due to fabrication variation, the response of each leg is slightly different, so the measurements were used to reconstruct the spring constant and the damping coefficient of each leg. For this study, only the first two resonant frequencies of the legs are considered in the model simulation. The relationship between vertical and lateral motion of the legs is also characterized for the four outer legs by measuring dynamic motion in the lateral direction with the LDV. A similar identification process is applied to the 30 mm prototype.

#### Coefficient of Friction and Restitution.

The contact properties of two types of materials (wood and metal) were tested as the ground for the robot prototypes. The coefficient of friction is determined for surfaces of different materials using at least five sliding tests per surface (Table 4). The coefficient of restitution is determined via measurement with LDV. The prototype is trapped by placing a large mass in front of robot to prevent motion in all other directions besides the z-direction (vertical). Sudden changes in z-direction velocity are recorded at least five separate times for each prototype leg, and all data are combined to calculate the average value of coefficients (Table 4).

### Model Results Versus Experimental Results.

To validate the dynamic model, the simulation results for leg and body motion based on the parameters in the “Model Validation” section are compared to the experimental results for both 80 mm and 30 mm prototypes under several forcing conditions. The 80 mm prototype is used primarily to validate the vertical motion of the legs and overall walking speed, and the 30 mm prototype is used primarily to validate the vertical motion of the body. The walking speed of each robot is characterized by actuating near resonance and recording the resulting motion with a camera. From the videos, the position of the robot is tracked in each frame, and the numerical derivative is computed to calculate the average forward velocity, with error range. These measurements are then compared to simulation results. Supplementary Video S1, which is available under the “Supplemental Materials” tab on the ASME Digital Collection, shows a representative sample of robot walking and subsequent velocity measurements for both prototypes. The validation study begins with the robot leg motion in the vertical direction, using the 80 mm prototype, where body motion is largely uniform. Robot body motion in vertical direction then shows some interesting phenomena in the 30 mm prototype, which is discussed next. Finally, the robot walking velocity predicted from the dynamic model is assessed using the 80 mm prototype, which is less effected by external wiring than the 30 mm robot. The comparison shows the influence from many design parameters toward the robot motion.

#### Leg Vertical Motion.

The time-domain contact responses of one leg of the 80 mm prototype under different input voltages are shown in Fig. 6, from experiments and simulation. The robot model accurately captures the general trends in contact duration, timing, and foot step amplitude in the vertical direction, for which high-fidelity measurement can be made directly with the LDV. The largest source of disagreement is the variation of motion amplitude at high-voltage input, which may be influenced by body vibration modes not included in the current model.

#### Robot Body Motion.

The modeling of the body's vertical motion is validated primarily using the 30 mm robot prototype. The 80 mm prototype is not well-suited for this since it has a larger mass, which leads to very small body motions. Characteristic time-domain responses of the vertical body motion of the 30 mm prototype are shown in Fig. 7. The measurement was performed with a square wave actuation signal at 1400 Hz with an amplitude of 30 V. Notably, as robot mass is reduced, the manner in which the robot's feet and ground interact is changed dramatically. Rather than producing foot impact at each actuation cycle (which occurred consistently for the 80 mm prototype), both simulation and experiments show that the body of the robot will stay in air for several leg cycles between foot impacts. Due to the extended time the body spends in air, the air damping term in Eq. (12) is considered for the 30 mm prototype. However, simulations for 80 mm prototype show that the air damping is negligible due to large prototype weight and limited vertical in-air motion. When the experimental and simulated results are compared, the amplitudes match for both the fast leg and slower body oscillations. To better quantify the in-air phase of the motion, a fast Fourier transform (FFT) is performed on the time-series data, as shown in Fig. 8. The first peaks of the FFT from the simulation and experiments are both located within the range of 50–100 Hz. In other words, the model successfully predicts the approximate range and frequency of impact events, even as robot scaling leads to more complex interactions between the robot and the ground. This phenomenon is similar to behavior previously observed in a magnetoelastic impact motor with similar modeling applied, which also showed bouncing of the rotor every few cycles of stator vibration [7].

#### Robot Velocity.

Locomotion experiments are designed to measure the relationship between actuation voltage, actuation frequency, ground condition, and average robot velocity. Due to significant influence of the power cord on the forward motion of the 30 mm prototype, forward velocity of the 80 mm prototype is primarily examined. However, it should be noted that simulations of the 30 mm robot yield a forward walking speed of 4 mm/s at 15 V actuation amplitude and 19 mm/s at 60 V, both of which lie inside the experimentally observed velocity ranges (2.2–4.1 mm/s at 15 V amplitude and 12–26 mm/s at 60 V).

The forward velocity of the 80 mm prototype was characterized using input frequencies primarily between 110 Hz and 140 Hz, because the first resonance modes of all six legs are located close to or in this range. The input voltage amplitude ranged from 15 to 45 V. Less than 15 V produced minimal movement from the robot, and 45 V was the upper limit of the experiment setup. Simulations using the dynamic model were also generated in these ranges. Comparison between simulations and experimental results shows that the direction of robot motion is the same. Error in simulation results is generated by ten simulations with randomly generated individual leg spring constants and motion amplitudes. The leg parameters are distributed according to the variability in the motion of experimentally observed legs, as measured during testing without ground contact. Since the range of leg parameters covers all the measured values of spring constants and motion amplitudes for all six legs, the error calculated is enough to predict the range of possible forward velocities. Error in experimental walking speed is measured from a 10–15 s robot walking video under each actuation condition. These error bars are intended to illustrate the potential variability in robot locomotion given the limitations of the accuracy with which nominal robot leg dynamics could be identified.

Figure 9 shows the simulated and measured velocities on a wood surface under different actuation frequencies, when the actuation voltage was fixed at 30 V. Figure 10 compares the behavior on the metal surface, showing representative trends with respect to voltage, when the actuation frequency was fixed at 140 Hz. The 140 Hz results are similar to the results around resonance (from 110 to 140 Hz). Predictions for robot velocity generally capture key trends in robot performance, though behavior is not as closely captured as was observed for vertical motion alone. Most importantly, velocity trends near the primary resonant peak are well-captured, but unmodeled dynamics take on more substantial roles at other frequencies. With respect to voltage, both experiments and simulation predict a local maximum for velocity. The existence of this local maximum is attributed to body motion starting to influence leg interaction, as shown in the previous leg contact measurements (Fig. 6). Body motion breaks up the timing between sets of legs in the tripod gait, but speed increases again as a function of voltage when individual leg displacements become sufficiently large to overcome this effect.

The importance of including the second resonance mode in the dynamics is also illustrated in Fig. 10. The trend in velocity with respect to actuation voltage follows a completely different trajectory than experiments when simulated using only the first resonance mode. The addition of the second resonance mode to the simulation resulted in more accurate dynamics estimation. Additionally, higher frequency modes (from third resonance) contribute a negligible portion of the overall motion for the robot size scales investigated, especially for displacement. If a hypothetical robot design did show sensitivity to higher modes, it should be noted that the dynamic model in this work is capable of including any further resonance.

The path of both prototypes is curved in experimental tests since the legs are not perfectly symmetrical in the fabricated prototypes. Simulation results shown in this paper assume uniform leg geometries but require material properties be measured experimentally for good agreement with walking experiments. Including leg variation has a very limited effect on forward and vertical velocity predictions (less than 10%), but can predict direction of curved walking paths (though with less accuracy).

These results do highlight some of the limitations of the proposed model. First, the discrepancy of the velocity in the frequency domain mainly occurs around 180 Hz, as shown in Fig. 9, where an additional peak appeared for a subset of the legs (Fig. 5). This peak is measured to be a compliant mode of the robot chassis, as opposed to vibration modes local to the individual legs. Conceptually, treating the entire robot as a compliant body would allow this mode to be captured and improve predictions of robot behavior. However, this increase in accuracy would come at a substantial increase in model complexity due to a larger number of modes to be accounted for at each leg and difficulty of analytically modeling chassis deformation. Also, the rigid chassis assumption is more appropriate for future implementation of payload, such as power and control units. Body resonance behavior changes with the presence of a payload, which should be calibrated case by case. Second, with respect to voltage, while the qualitative behavior of local maxima is successfully captured and velocities are mostly within the margin of error, the predicted locally optimal voltage is significantly different. This could be due to neglecting nonlinearities in the piezoelectric actuator response or some adhesion effects at the robot foot.

In addition to actuation voltage and frequency, the motion amplitude ratio between vertical and axial motion was examined, as shown in Fig. 12. The ratio of y- and z-direction motion is determined as the maximum free motion amplitude in z-direction over y-direction. Given a fixed vertical motion, an increase in axial motion is correlated with an increase in the robot's forward velocity, as expected, but with a more pronounced local optimum. These results indicate the importance of characterizing the correct ratio between two directional motions. The influence of variations in spring constants and damping coefficients of each leg are shown as error bars in Figs. 912. Higher actuation voltage generally produces a large error range since the absolute value of the motion amplitudes is larger with the same error range. Also, the error in actuation frequency away from resonance is larger than the one around resonance. Moreover, the robot is able to move in the reverse direction as indicated by simulation. This phenomenon was also observed in experiments, but not quantitatively measured.

Similarly, the effect of foot–terrain interaction parameters is shown in Fig. 13. The robot's velocity drops substantially if the coefficient of friction drops under 0.1, which indicates the existence of slippage between the robot foot and ground. Above this threshold, the exact coefficient of friction has little influence, as the stationary friction force between the feet and ground is nearly always less than the maximum friction force that could be provided. Meanwhile, low coefficients of restitution are predicted to lead to higher robot velocities, because lower vertical velocities after contact will cause the foot to stay in contact with ground longer. Longer ground contact allows friction to actuate the robot forward, provided that the in-plane and out-of-plane motion are coupled in the axes assumed by this work. Low coefficients of restitution are also related to plastic contact between the foot and ground. Regarding the coefficient of restitution and friction, both types of materials we used as ground condition share close value.

## Discussion and Conclusion

A model is presented for dynamics of small, elastic walking robots. The model includes 3DOF rigid body motion and 2DOF flexible leg motion with multiple modes to capture trends in robot leg motion and velocity. The interaction between the robot and the ground is modeled using simple coefficients of friction and restitution. To validate this model, two different centimeter-scale robot prototypes are fabricated, both of which can be represented by the model. The prototypes are actuated by piezoelectric ceramic strips bonded to a 3D-printed body and are experimentally characterized in terms of leg and body dynamics as well as walking speed. The velocities of the robot prototypes under varying actuation voltages and frequencies near the legs' resonant frequencies are measured to be around 3 mm/s for both prototypes, which are in the same range as those simulated by the dynamic model after robot parameters are identified experimentally. The payload effect is also considered and validated with this model.

Other experimentally observed features that are captured well by the model include: (1) vertical leg motion trajectories and ground interaction trends; (2) bouncing effects, where large vertical body motion has multiple significant frequency components at frequencies lower than the leg actuation frequency, as exhibited in the 30 mm prototype; (3) velocity trends with changing voltage, when a sufficient number of vibration modes are accounted for in the leg model; and (4) the existence of a local maximum velocity with respect to voltage at some frequencies, though the maximum occurs at slightly higher voltages in simulation than in experiments. Simulated robot velocity trends with respect to the robot legs' vertical–axial motion ratio, CoF and CoR are also reported. Limitations of the model include omitting the possibility for body deformation to simplify the model and the need for a significant amount of leg and terrain model information to accurately predict the robot behavior.

As a structure for future microrobot development, this work provides: (1) a generalized model for multiple-legged systems that includes both rigid body and compliant actuation motion; (2) accommodation of leg motion with multiple modes and DOF while body motion has multi-DOF; and (3) good performance with high-frequency, lightly damped actuation. These features are anticipated to be common in future millimeter-scale robots based on microfabrication processes, when operated in dynamic walking or running gaits. Other phenomena that may be added in future work for even smaller robots include effects such as squeeze-film damping, adhesion or stiction forces at robot feet, and other small-scale forces.

## Acknowledgment

This work was supported by NSF awards CMMI 1435222, IIS 1208233, and CMMI 0954422. The authors thank Mr. Jongsoo Choi for the fabrication of the millimeter-scale millipede robot shown as motivation and for additional assistance with robot prototyping and testing.

## References

1.
Kim
,
S.
,
Clark
,
J. E.
, and
Cutkosky
,
M. R.
,
2006
, “
iSprawl: Design and Tuning for High-Speed Autonomous Open-Loop Running
,”
Int. J. Rob. Res.
,
25
(
9
), pp.
903
912
.
2.
Cham
,
J. G.
,
Bailey
,
S. A.
,
Clark
,
J. E.
,
Full
,
R. J.
, and
Cutkosky
,
M. R.
,
2002
, “
Fast and Robust: Hexapedal Robots Via Shape Deposition Manufacturing
,”
Int. J. Rob. Res.
,
21
(
10–11
), pp.
869
882
.
3.
Uemura
,
M.
, and
Kawamura
,
S.
,
2009
, “
Resonance-Based Motion Control Method for Multi-Joint Robot Through Combining Stiffness Adaptation and Iterative Learning Control
,”
IEEE International Conference on Robotics and Automation
(
ICRA
), Kobe, Japan, May 12–17, pp.
1543
1548
.
4.
Ruderman
,
M.
,
Hoffmann
,
F.
, and
Bertram
,
T.
,
2009
, “
Modeling and Identification of Elastic Robot Joints With Hysteresis and Backlash
,”
IEEE Trans. Ind. Electron.
,
56
(
10
), pp.
3840
3847
.
5.
Hoffman
,
K. L.
, and
Wood
,
R. J.
,
2011
, “
Passive Undulatory Gaits Enhance Walking in a Myriapod Millirobot
,”
IEEE/RSJ International Conference on Intelligent Robots and Systems
(
IROS
), Taipei, Taiwan, Oct. 18–22, pp.
1479
1486
.
6.
Hoffman
,
K. L.
, and
Wood
,
R. J.
,
2011
, “
Myriapod-Like Ambulation of a Segmented Microrobot
,”
Auton. Rob.
,
31
(
1
), pp.
103
114
.
7.
Qu
,
J.
,
Tang
,
J.
,
Gianchandani
,
Y. B.
,
Oldham
,
K. R.
, and
Green
,
S. R.
,
2015
, “
Dynamic Modeling of a Bidirectional Magnetoelastic Rotary Micro-Motor
,”
Sens. Actuators A
,
223
, pp.
49
60
.
8.
Yasuda
,
T.
,
Shimoyama
,
I.
, and
Miura
,
H.
,
1994
, “
Microrobot Locomotion in a Mechanical Vibration Field
,”
,
9
(
2
), pp.
165
176
.
9.
,
P. E.
, and
Bright
,
V. M.
,
2000
, “
Prototype Microrobots for Micro-Positioning and Micro-Unmanned Vehicles
,”
Sens. Actuators A
,
80
(
2
), pp.
132
137
.
10.
Murthy
,
R.
,
Das
,
A. N.
,
Popa
,
D. O.
, and
Stephanou
,
H. E.
,
2011
, “
ARRIpede: An Assembled Die-Scale Microcrawler
,”
,
25
(
8
), pp.
965
990
.
11.
Shin
,
M.
,
Choi
,
J.
,
Rudy
,
R. Q.
,
Kao
,
C.
,
Pulskamp
,
J. S.
,
Polcawich
,
R. G.
, and
Oldham
,
K. R.
,
2014
, “
Micro-Robotic Actuation Units Based on Thin-Film Piezoelectric and High-Aspect Ratio Polymer Structures
,”
ASME
Paper No. DETC2014-35145.
12.
Rhee
,
C.-H.
,
Pulskamp
,
J. S.
,
Polcawich
,
R. G.
, and
Oldham
,
K. R.
,
2012
, “
Multi-Degree-of-Freedom Thin-Film PZT-Actuated Microrobotic Leg
,”
J. Microelectromech. Syst.
,
21
(
6
), pp.
1492
1503
.
13.
Ryou
,
J. H.
, and
Oldham
,
K. R.
,
2014
, “
Dynamic Characterization of Contact Interactions of Micro-Robotic Leg Structures
,”
Smart Mater. Struct.
,
23
(
5
), p.
055014
.
14.
Firebaugh
,
S.
,
Piepmeier
,
J.
,
Leckie
,
E.
, and
Burkhardt
,
J.
,
2011
, “
Jitterbot: A Mobile Millirobot Using Vibration Actuation
,”
Micromachines
,
2
(
2
), pp.
295
305
.
15.
Casarez
,
C. S.
,
Ryou
,
J. H.
, and
Oldham
,
K. R.
,
2013
, “
Dimensional Analysis of Dynamic MEMS Micro-Robotic Walking Subject to Orthogonal Actuation and Small-Scale Forces
,”
ASME
Paper No. DETC2013-13473.
16.
Qu
,
J.
, and
Oldham
,
K. R.
,
2016
, “
Modeling Legged Micro-Robot Locomotion Based on Contact Dynamics and Vibration in Multiple Modes and Axes
,”
ASME
Paper No. IDETC2016-59621.
17.
Saranli
,
U.
,
Buehler
,
M.
, and
Koditschek
,
D. E.
,
2001
, “
RHex: A Simple and Highly Mobile Hexapod Robot
,”
Int. J. Rob. Res.
,
20
(
7
), pp.
616
631
.
18.
Birkmeyer
,
P.
,
Peterson
,
K.
, and
Fearing
,
R. S.
,
2009
, “
DASH: A Dynamic 16g Hexapedal Robot
,”
IEEE/RSJ International Conference on Intelligent Robots and Systems
(
IROS
), St. Louis, MO, Oct. 10–15, pp.
2683
2689
.
19.
Lambrecht
,
B. G.
,
Horchler
,
A. D.
, and
Quinn
,
R. D.
,
2005
, “
A Small, Insect-Inspired Robot That Runs and Jumps
,”
IEEE International Conference on Robotics and Automation
(
ICRA
), Barcelona, Spain, Apr. 18–22, pp.
1240
1245
.
20.
Morrey
,
J. M.
,
Lambrecht
,
B.
,
Horchler
,
A. D.
,
Ritzmann
,
R. E.
, and
Quinn
,
R. D.
,
2003
, “
Highly Mobile and Robust Small Quadruped Robots
,”
IEEE/RSJ International Conference on Intelligent Robots and Systems
(
IROS
), Las Vegas, NV, Oct. 27–31, pp.
82
87
.
21.
Hoover
,
A. M.
,
Steltz
,
E.
, and
Fearing
,
R. S.
,
2008
, “
RoACH: An Autonomous 2.4 g Crawling Hexapod Robot
,”
IEEE/RSJ International Conference on Intelligent Robots and Systems
(
IROS
), Nice, France, Sept. 22–26, pp.
26
33
.
22.
Brufau-Penella
,
J.
,
Sánchez-Martín
,
J.
, and
Puig-Vidal
,
M.
,
2006
, “
Piezoelectric Polymer Model Validation Applied to mm Size Micro-Robot I-SWARM (Intelligent Swarm)
,”
Proc. SPIE
,
6166
, p.
61660Q
.
23.
Woern
,
H.
,
Szymanski
,
M.
, and
Seyfried
,
J.
,
2006
, “
The I-SWARM Project
,”
The 15th IEEE International Symposium on Robot and Human Interactive Communication
(
ROMAN
), Hatfield, UK, Sept. 6–8, pp.
492
496
.
24.
Ebefors
,
T.
,
Mattsson
,
J. U.
,
Kälvesten
,
E.
, and
Stemme
,
G.
,
1999
, “
A Walking Silicon Micro-Robot
,” The 10th International Conference on Solid-State Sensors and Actuators (
Transducers '99
), Sendai, Japan, June 7–10, pp.
1202
1205
.
25.
DeVoe
,
D. L.
, and
Pisano
,
A. P.
,
1997
, “
Modeling and Optimal Design of Piezoelectric Cantilever Microactuators
,”
J. Microelectromech. Syst.
,
6
(
3
), pp.
266
270
.
26.
Smits
,
J. G.
, and
Choi
,
W.-S.
,
1991
, “
The Constituent Equations of Piezoelectric Heterogeneous Bimorphs
,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Control
,
38
(
3
), pp.
256
270
.
27.
Rao
,
S. S.
, and
Yap
,
F. F.
,
1995
,
Mechanical Vibrations
,
,
New York
.
28.
Gilardi
,
G.
, and
Sharf
,
I.
,
2002
, “
Literature Survey of Contact Dynamics Modelling
,”
Mech. Mach. Theory
,
37
(
10
), pp.
1213
1239
.