Optimal Gaits for Drag-dominated Swimmers with Passive Elastic Joints
OOptimal Gaits for Drag-dominated Swimmers with Passive Elastic Joints
Suresh Ramasamy ∗ and Ross L. Hatton † Collaborative Robotics and Intelligent Systems (CoRIS) Institute,Oregon State University, Corvallis, OR, USA (Dated: October 5, 2020)In this paper, we identify optimal swimming strategies for drag-dominated swimmers with apassive elastic joint. We use resistive force theory (RFT) to obtain the dynamics of the system.We then use frequency domain analysis to relate the motion of the passive joint to the motion ofthe actuated joint. We couple this analysis with elements of the geometric framework introducedin our previous work aimed at identifying useful gaits for systems in drag dominated environments,to identify speed-maximizing and efficiency-maximizing gaits for drag-dominated swimmers with apassive elastic joint.
I. INTRODUCTION
A common strategy for locomotion among animals androbots is to couple cyclical shape changes (gaits) to aninteraction with the environment. A long term researchfocus of the geometric mechanics community has beenfinding geometric principles that describe what makes agait effective. Many geometric tools resulting from thisline of research have only been developed for systemswith a fully actuated shape space. Passive elastic joints,such as flexible fins and tails, however, play an importantrole in locomotion in many systems [1]. In this paper,we expand on the geometric framework that was devel-oped in [2, 3] to identify optimal gaits for fully actuatedswimmers in low Reynolds number fluids, extending it toswimmers with passive elastic shape elements. We thenuse this extended framework to identify optimal gaits forswimmers in drag dominated environments with passiveelastic shape elements.Our analysis draws on a rich history of work in the ge-ometric mechanics community aimed at using conceptsfrom differential geometry to understand how systemslocomote in a low Reynolds number fluid, what the mainfactors that determine the efficiency of gaits are, andwhat optimal gaits look like for these systems. In par-ticular, we know that the efficiency of a gait (speed ata given power level or power required at a given speed)depends on its path, period, and pacing in the system’sshape space:1. The net displacement per cycle corresponds to theamount of “constraint curvature” the gait encom-passes. Gaits that maximize displacement per cy-cle thus enclose sign-definite regions of the system’sconstraint curvature functions (CCFs), which canbe obtained from the dynamics of the system [4, 5]. ∗ [email protected] † [email protected]
2. At a given level of average power consumption, theminimum time taken to execute one cycle of thegait corresponds to the path length of the gait tra-jectory under a Riemmanian metric on the shapespace, meaning that shorter gaits can be executedmore frequently or at a lower power level thanlonger gaits [6].3. Gaits that progress through their gait kinematicsat a steady pace have a lower average power costfor a given frequency than those that “surge” and“dwell” at points in the cycle, or equivalently, canbe executed at a higher frequency for a given aver-age power cost [7].The process of finding efficient gaits thus involves strik-ing a balance between maximizing the enclosed constraintcurvature and minimizing the metric-weighted perimeterof the gait and then finding a steady-pace parametriza-tion of the resulting curve. In [2, 3, 8], we built a vari-ational framework based on this geometric insight thatidentifies optimal gaits for fully-actuated drag-dominatedkinematic systems.In this paper, we extend this framework to identify op-timal gaits for swimmers with passive elastic joints anddemonstrate the framework on systems with a single ac-tive and a single passive joint. The ways in which thedynamics of swimmers with passive elastic joints differfrom the dynamics of fully actuated swimmers are:1. Because of the coupling between the actuated andpassive joints, the passive swimmers can executeonly a subset of the gait kinematics the fully actu-ated swimmers can execute.2. This coupling between the actuated and passivejoints also endows each gait achievable by the pas-sive swimmer with a unique pacing. Hence, thepacing cost cannot be minimized separately fromthe kinematic cost.As illustrated in Fig. 1, we address the two problemsintroduced by the presence of passive elastic joints: wefirst use frequency domain analysis to analytically ap-proximate the motion of the unactuated joint in response a r X i v : . [ phy s i c s . b i o - ph ] S e p FrequencyFrequency P h a s e d i ff e r e n c e (a) Optimal gait for the fully actuated Purcell Swimmer (b) Frequency response of the passive joint to oscillations of the active joint (c) Optimal gait for the passive Purcell swimmer (rad) ( r a d ) (rad) ( r a d ) (rad/s)(rad/s)(deg) FIG. 1: Effect of a passive elastic joint on the shape and pacing of the (speed-maximizing) optimal gait for a Purcell swimmer. (a) The(speed and efficiency maximizing) optimal gait for the fully actuated Purcell swimmer overlaid on the curvature of the system dynamics.(b) Illustration of the Bode plot of the response of the passive elastic joint to oscillations of the active joint. This response dictates thelocus of achievable gaits in the shape space for the passive swimmers. (c) The (speed maximizing) optimal gait for the passive Purcellswimmer. Whereas the optimal gait for the fully actuated Purcell swimmer is plotted with a line of uniform thickness indicating constantpower dissipation throughout the cycle, the optimal gait for the passive Purcell swimmer in (c) is plotted using a line of varyingthickness, with thickness at a point corresponding to the magnitude of power required at that point of the gait. to the motion of the actuated joint. We then combine thisfrequency-space analysis with elements of the geomet-ric framework introduced in [8] to construct a gradient-descent algorithm that identifies optimal gaits and thepacing associated with these gaits for passive swimmers.The optimal gaits for passive swimmers maximize theCCF integral relative to perimeter and pacing costs, sub-ject to amplitude and phase constraints of a first ordersystem.Drag-dominated swimmers with passive elements werepreviously studied in [9–13]. Of these works, [9, 10]and [11] are most relevant to this paper, as they dis-cuss the motion of swimmers with a harmonically-drivenactive joint and a passive joint. The analysis in [9] isparticularly relevant, where using perturbation expan-sion, explicit expressions for leading order solutions werederived for harmonic input oscillations, and the optimalswimmer geometry was obtained for the Purcell swim-mer. In this paper, unlike in [9], we do not restrict ourinput to simple harmonic oscillations, and use a higherorder representation of the system dynamics.The rest of the paper is organized as follows: In § II,we review the geometric locomotion model and derivethe dynamics of the drag-dominated swimmers. In § III,we present our frequency domain analysis to relate themotion of the passive joint to the motion of the activejoint. In § IV, we recall a few key elements of the gra-dient descent algorithm introduced in [8], and combinethem with the frequency domain analysis of § III, to setup the stroke optimization problem. In § V, we describethe process of finding the speed-maximizing gaits for pas- sive swimmers. In § VI, we describe the process of findingthe efficiency-maximizing gaits for the passive swimmers.
II. SYSTEM DYNAMICS
In this section, we review the dynamics of swimmersin a low Reynolds number fluid. We focus here on thegeometric structure of the system dynamics, and not onthe explicit expressions. A more detailed discussion ofhow the dynamics are obtained is presented in [4] and [6].Some of the the materials presented in this section alsoappear in [4] and [6].
A. System geometries
The Purcell swimmer illustrated in Fig. 2(a) is a clas-sic example three-link system with minimal complexitythat can swim in a low Reynolds number fluid (drag-dominated environment), and was introduced in [14].The T-link swimmer illustrated in Fig. 2(b) is a modifi-cation of the Purcell swimmer where one of the peripherallinks is attached to the center link at its midpoint. TheT-link swimmer was introduced in [10] as a simplifiedmodel to study the swimming dynamics of
Schistosomamansoni . S.mansoni causes schistosomiasis, a diseasecomparable to malaria in socio-economic impact [10].In the passive Purcell swimmer, a torsional spring isattached at the second joint as shown in Fig. 2(a) suchthat when the joint angle α is zero, the spring is at itsequilibrium. In the passive T-link swimmer, the torsionalspring is attached between the middle link and the pe-ripheral link attached at its midpoint as shown in 2(b)such that when the joint angle α is zero, the spring isat its equilibrium. B. Geometric Locomotion Model
When analyzing a locomoting system, it is convenientto separate its configuration space Q (i.e. the space of itsgeneralized coordinates q ) into a position space G anda shape space R , such that the position g ∈ G locatesthe system in the world, and the shape r ∈ R gives therelative arrangement of the particles that compose it [15].The locomotion model we employ in this paper was de-veloped for systems in kinematic regimes where no glid-ing can occur i.e, where zero shape velocity results in zeroposition space velocity [16]. In this model, there exists alinear relationship at each shape between changes in thesystem’s shape and changes in its position, ◦ g = − A ( r ) ˙ r, (1)in which ◦ g = g − ˙ g is the body velocity of the system(i.e., ˙ g expressed in the system’s local coordinates), andthe local connection A linearly maps joint velocities tothe body velocity they produce by pushing the systemagainst its environment. Each row of − A can be re-garded as a body-coordinates local derivative of one po-sition component with respect to the system shape. Ifwe plot the rows of − A as arrow fields, as in Fig. 3, thismeans that moving in the direction of the arrows movesthe system positively in the corresponding body direc-tion, and moving perpendicular to the arrows results inno motion in that direction [5, 17].Several efforts in the geometric mechanics community(including our own), have aimed to use the structureof the systems’ Lie brackets (a measure of how “non-canceling” the system dynamics are over cyclic inputs)to understand the structure of the optimal solutions tothe system equations of motion [4, 6, 18–23]. The coreprinciple in these works is that because the net displace-ment g φ over a gait cycle φ is the line integral of (1) along φ , the displacement can be approximated [24] by an inte-gral of the curvature D ( − A ) of the local connection (itstotal Lie bracket, which is a measure of how much theconnection changes, over a surface φ a bounded by thecycle) [5]: g φ = ‰ φ − g A ( r ) (2) ≈ exp ¨ φ a − dA + (cid:88) i In this paper, as in [4, 6], we generate the dynamicsfor our example systems from a resistive force model, inwhich each element of the body is subject to normal and (a) Connection vectorfield (b) Constraint CurvatureFunction FIG. 3: The connection vector fields in the x-direction for thePurcell swimmer and T-link swimmer are shown in (a) and thecorresponding constraint curvature functions are shown in (b). [4]. tangential drag forces proportional to their velocities inthose (local) directions [26]. The normal drag coefficientis larger than the tangential component (here, by a factorof 2 : 1), corresponding to the general principle that it isharder to move a slender object in a fluid or on a surfacecrosswise than it is to move it along its length [27].We combine this linear drag condition with a quasi-static equilibrium condition that the net drag force andmoment on the system is zero at all times (treating thesystem as heavily overdamped, with acceleration forcesmuch smaller than drag forces). Because the drag forcesare not isotropic, the quasi-static condition does not pre-vent the system from moving, and the system can usethe angle-of-attack of its body surfaces to generate netmotion.The quasistatic condition imposes a Pfaffian con-straint [28] on the system’s generalized velocity, whichstates that combinations of body and shape velocitiesthat are consistent with the system physics must lie inthe null space of a linear map from the system velocitiesto the forces and moments acting on the system’s bodyframe, F bx F by F bθ = = ω ( r ) (cid:20) ◦ g ˙ r (cid:21) . (6)The matrix ω that maps the velocities to the net forces onthe body frame is a pullback of the drag matrices of theindividual links via the system kinematics and dependsonly on the shape r . By separating ω into two sub-blocks, ω = [ ω × g , ω × nr ],it is straightforward to rearrange (6) into ◦ g = − ( ω − g ω r ) ˙ r, (7)revealing the local connection as A = ω − g ω r . The ex-pressions for the dynamics are unwieldy (running to sev-eral pages of trigonometric terms in even the simplestcases) so we do not write them out in full here. See [4]for a more detailed treatment of (6)–(7) in the case ofthe three-link swimmer).Using the local connection A and the system’s internalkinematics, we can obtain the Jacobian from shape veloc-ity to the local velocity of each section of the body J ( r, (cid:96) ),where (cid:96) is the location of the section on the body. Wecan use this Jacobian to calculate a Riemannian metric M over the shape space that encodes the cost of effectinga shape change as M ( r ) = ˆ body (cid:16) J T ( r, (cid:96) ) (cid:17) C (cid:16) J ( r, (cid:96) ) (cid:17) d(cid:96), (8)where C is the matrix of drag coefficients, which actsas a local metric for the motion of each element of thebody. For the systems considered in the paper, we takethis matrix as, C = (9)indicating that for any infinitesimal element of a link,the resistance to lateral motion is twice the resistanceto longitudinal motion and that the moment-drag on alink is generated by the translation of the longitudnal andlateral forces acting in the local frame of the infinitesimalelement to the body frame of the system [29].As discussed in our previous work [4, 5], the metric M encodes a quadratic relationship between the shape ve-locities and power dissipated into the surroundings, givenby P = ˙ r T M ( r ) ˙ r, (10)as well as the mapping from joint velocities to torques onthe joints, τ = M ( r ) ˙ r. (11)Details of the calculations to generate the local connec-tion A and the Riemmanian metric M for our examplesystems are provided in [4] and [6].In drag-dominated environments, a common measureof the cost of any motion executed by a swimmer is theenergy dissipated into the surrounding fluid while exe-cuting the motion, E = ˆ T P ( t ) dt (12)Substituting the expression for power dissipatedfrom (10) provides a measure of the cost in terms of themetric and shape trajectory, E = ˆ T P ( t ) dt = ˆ T ˙ r T M ( r ) ˙ rdt. (13)This cost depends on the geometry, time period and pac-ing of a given gait. As discussed in [7, 30, 31], for fullyactuated swimmers this cost can be broken down into acombination of a pacing-invariant cost that measures thepathlength s of the trajectory through the shape space(weighted by the shapespace metric M ), and a pacingcost σ that measures the deviation from optimal pac-ing. Finding a gait that minimizes the energy dissipatedinto the surrounding fluid E is thus equivalent to find-ing a gait the minimizes the metric-weighted pathlength s , and executing said gait at a constant-power pacing tominimize σ .The pathlength s of a curve r ( t ) under a metric M is s = ˛ ds = ˛ (cid:113) dr T M ( r ) dr. (14)Changing the variable of integration from shape r to time t , we can relate the pathlength to the square root of powerexpended: s = ˆ T (cid:113) ˙ r T M ( r ) ˙ rdt = ˆ T (cid:112) P ( t ) dt (15)Because moving with constant power is the least-costlypacing with which to execute a motion under viscousdrag [7], we can further simplify (15) as s = (cid:112) P avg T, (16)where P avg is the average power utilized while executingthe motion. This pathlength provides a geometric costfor the best-case execution of the kinematics in a gaitcycle.The additional cost for a non-optimal pacing can berepresented by squaring the difference between the aver-age and instantaneous rates at which the gait is beingfollowed (measured as s per time), and then integratingover the time during which the gait is being executed, σ = ˆ τ total (cid:18) s total τ total − ddτ ( s ( τ )) (cid:12)(cid:12)(cid:12) τ = t (cid:19) dt, (17)where τ total is the time period of the gait, s total is thelength of the gait under the metric M , and s is distancetraveled along the gait as a function of time correspond-ing to the given pacing. If the gait is proceeding at con-stant power, s total τ total is equal to the rate at which s changeswith time, so σ measures the extent to which the pacinglags and leads the optimal pacing.Any pacing other than constant-power will make thetrajectory take longer for a given average power (or in-crease the average power required to complete the motionin a fixed time). III. FREQUENCY DOMAIN ANALYSIS The key difference in the dynamics of a swimmer with apassive joint when compared to a fully actuated swimmeris the coupling of the motion of the actuated and unac-tuated joint. In this section, we explore this differencefurther and present a way of accurately approximatingthe motion of the unactuated joint from the motion ofthe actuated joint using frequency domain analysis. Themethod of linearizing the passive dynamics to obtain ap-proximate limit cycles presented in this section is in thesame vein as the limit cycle analysis presented in [13],where a two-link system with static separation betweencenters of mass and buoyancy was studied. A. Dynamics of the passive elastic joint As discussed in [4, 5], the mapping between joint ve-locities and torques on the joint in the fully actuatedswimmers is encoded in the metric calculated in (8), τ = M ( r ) ˙ r. (18)Because the systems considered in this paper have onlyone active and one passive joint, this relationship be-comes (cid:20) τ τ (cid:21) = M ( α , α ) (cid:20) ˙ α ˙ α (cid:21) . (19)In the case of the swimmers with an elastic joint, becausethe actuation in the second joint is replaced by an elasticelement with stiffness k , the torque τ is always equal to − kr , i.e., (cid:20) τ − kα (cid:21) = M ( α , α ) (cid:20) ˙ α ˙ α (cid:21) . (20)The first equation in this system of equations, τ = M ˙ α + M ˙ α , (21)thus relates the torque in the actuated joint to the mo-tion of the joints and can be used to calculate the torquerequired to effect any feasible motion. The second equa-tion in this system of equations, − kα − M ˙ α = M ˙ α (22)or equivalently, − k M α − M M ˙ α = ˙ α (23)encodes the dynamics of the passive elastic joint in termsof the active joint, and thus defines the space of feasiblemotions.Since the joint α is assumed to be the actuated joint,we have full control of ˙ α . The value of M depends on α and α . Its dependence on α and α conveys how FIG. 4: Comparison of the exact stroke limit cycles with theshape predicted by frequency domain analysis. Each black solidline represents the motion of the swimmer when the input to theactuated joint is a sinusoidal wave. We can see the motionconverges to a limit cycle. The red dashed line represent theshape of the limit cycles predicted by frequency domain analysispresented in § III. The cartoons in the background show how thePurcell swimmer looks like at different points of the shape space. the shape of the robot affects the effort required to movethe joints. If we consider gaits that are relatively smalloscillations of shape, we can approximate the value of M to be constant throughout the gait. This assumptionnecessarily introduces errors in our prediction of the mo-tion of the passive joint when the amplitude of input tothe active joint is large.For both the Purcell and T-link swimmers with passiveelastic joints, assuming the value of M to be constantdoes not introduce significant errors for gaits of ampli-tude up to 1.5 radians: In Fig. 4, we illustrate the dis-tortions caused in the shape of the limit cycles when weassume M to be constant throughout the shape space.Each solid black line represents the motion of the fullswimmer model when the input to the actuated joint is asinusoidal wave, and the system starts with both anglesat zero. There is a transient term that dominates beforethe system reaches the limit cycle. The red dashed linesrepresent the shape of the limit cycles predicted when weassume M is constant (we describe the calculation of the K eq C eq FIG. 5: Spring damper system whose dynamics are equivalent tothe passive dynamics of the Purcell swimmer with a passiveelastic joint. limit cycle in the next subsection). B. Transfer function analysis In this subsection, our goal is to obtain an analyti-cal approximation of the response of the passive jointto input oscillations of the actively controlled joint. Weassume M to be constant throughout the gait, whichmakes (23)— the equation that describes the dynamicsof the passive elastic joint— a linear first order differen-tial equation and thus well suited to frequency domainanalysis.We use the Laplace transform on (23) to obtain thedynamics of the passive elastic joint in the frequency do-main as K eq L ( α ( t )) + C eq s L ( α ( t )) = s L ( α ( t )) (24)where K eq = − k M and C eq = −M M . These substi-tutions reveal that the dynamics of the passive joint areequivalent to those of a massless particle attached to afixed base through a spring and being driven through adamper by a position trajectory α as shown in Fig. 5.Note that the damping coefficient is completely depen-dent on the physics of the interaction between the swim-mer and the fluid.We can rewrite (24) as a transfer function relation be-tween the active joint α and passive elastic joint α L ( α ( t )) = H ( s ) L ( α ( t )) (25)where H ( s ) = s ( C eq s + K eq ) is the transfer function thatencodes the response of the passive elastic joint to oscial-lations of the active joint.Equation (25) tells us how inputs to the controlledjoint, α are mapped to the response of the passive jointin the frequency domain. In order to find the responseof the passive joint to a sinusoidal oscillation of the ac-tuated joint, we let α ( t ) = sin( wt ), then using (25), we Change in amplitude of input Change in frequency of input(b) Amplitude=1, Frequency=0.5(a) Amplitude=1.2, Frequency=0.5 (c) Amplitude=1, Frequency=0.25 FIG. 6: Changes in the shape of a gait resulting from changes to amplitude and frequency of input oscillation. (a) The gait resulting froma change in amplitude of the nominal input (black) is a scaled version of the gait corresponding to the nominal input actuation (red). (b)The gait resulting from nominal input actuation. (c) The gait resulting from a change in the frequency of the nominal input (black). Achange in the frequency of the nominal input leads to a change in both the amplitude and phase of the response of the passive joint. obtain L ( α ( t )) = H ( s ) (cid:122) (cid:125)(cid:124) (cid:123) s ( C eq s + K eq ) L ( α ( t )) (cid:122) (cid:125)(cid:124) (cid:123) ω ( s + ω ) (26)= A ( C eq s + K eq ) (cid:124) (cid:123)(cid:122) (cid:125) Transient term + A s + A ( s + ω ) (cid:124) (cid:123)(cid:122) (cid:125) Phase shifted sine wave (27)where A , A and A can be obtained by equating thetwo expressions for L ( α ( t )) in (26) and (27) as sω = A ( s + ω ) + ( A s + A )( C eq s + K eq ) , (28)and equating the coefficients of powers of s on each sideto extract a system of three equations, C eq K eq C eq ω K eq A A A = ω , (29)which can be easily solved to obtain A , A and A as A A A = C eq K eq C eq ω K eq − ω . (30)If our actuated joint follows a trajectory given by α ( t ) = B sin( ωt ), after reaching the steady state pe-riodic orbit, the unactuated joint, therefore, follows thetrajectory α ( t ) = B A cos( ωt ) + B A sin( ωt ) . (31)If our actuated joint follows a trajectory given by α ( t ) = B sin( ωt ) + B cos( ωt ), after reaching the steady stateperiodic orbit, the unactuated joint follows the trajectory α ( t ) = B ( A cos( ωt ) + A sin( ωt ))+ B ( − A sin( ωt ) + A cos( ωt )) . (32) Note that the value of A and A depend on the valueof ω , the frequency of the sinusoidal input to the actuatedjoint. The image of the gait in the shape space is thuscoupled to the pacing of the input to the actuated joint.This coupling between the image of the gait in theshape space and the input to the actuated joint is illus-trated in Fig. 6. Fig. 6(b) shows the shape of the gait re-sulting from a sinusoidal oscillation of the actuated jointof amplitude 1 and frequency 0.5. Fig. 6(a) shows theeffect of increasing the amplitude of oscillation on theshape of the gait. As we have linearized the dynam-ics of the passive joint, a change in amplitude withouta change in frequency produces a scaled version of theoriginal gait. Fig. 6(c) shows the effect of decreasing thefrequency of oscillation on the shape of the gait, with aslower oscillation leading to a weaker response (in termsof magnitude) from the passive joint.The Bode plot of the transfer function H in (25), pre-sented in Fig. 7, illustrates the specific nature of the rela-tionship. We can see that if the frequency of input to thecontrolled joint is very high, the response of the passiveelastic joint is phase shifted relative to the actuated jointby almost 180 degrees, resulting in a gait with a smallarea. At very low frequencies, the magnitude of the re-sponse of the passive elastic joint is small and results ina gait with a small area. At mid-range frequencies, how-ever, the response of the passive elastic joint has a largeramplitude and a phase shift that results in gaits withlarger areas. IV. FINDING OPTIMAL GAITS Optimal gait design has a long history of research inthe physics, mathematics, and engineering communities,as part of the broader field of optimal control [32, 33].Notable contributions to finding optimal gaits for swim-mers in a drag-dominated environment include those ofPurcell, who introduced the three-link swimmer as a min-imal template for understanding locomotion, a series of FIG. 7: Bode plot of the transfer function relating the output ofthe passive joint to the input of the actuated joint in the Purcellswimmer with a passive elastic tail. The inset figures show howthe periodic orbits corresponding to gain and phase at certainfrequencies look like in the shape space. We can see thatactuation at very low frequencies leads to gaits that enclose verylittle surface area due to the amplitude of the passive joint beinglow, while at very high frequencies the surface area enclosed bythe gait is low due to the almost 180 degree phase shift betweenthe oscillations of the actuated and passive joints. works [34–38] aimed at numerically optimizing the strokepattern, and the observation in [7] that the optimal pac-ing for the gait keeps the power dissipation constant overthe cycle.In drag-dominated environments, a common measureof the cost of any motion executed by a swimmer is thepower dissipated into the surrounding fluid while execut-ing the motion. A natural choice of definition for effi-ciency for these systems is thus η = g φ E , (33)where g φ is the displacement produced in the x directionover a gait cycle φ and E as defined in (12) is the totalenergy dissipated into the surrounding while executingthe gait φ .Note that this definition of efficiency is the inverse ofthe mechanical cost of transport used, e.g., in [30]. Themechanical cost of transport is a widely used efficiencymetric in the robotics community, especially while study-ing legged locomotion. As explained in Appendix A, gaits which optimize our criterion also optimize Lighthill’s ef-ficiency [39], which compares the power dissipated whileexecuting the gait to the power dissipated in rigidlytranslating the swimmer through the fluid; this measurehas an advantage over Lighthill in that it allows for ef-fective comparison between systems with different mor-phologies and does not require designating a referenceshape for rigidly dragging the swimmer.As noted in § II C, under optimal pacing (constant-power pacing), the pathlength s defined in (14) and thetotal power dissipated E are equivalent cost criteria. Asdetailed in [8], for fully actuated swimmers this insightlets us restrict our optimization to constant power tra-jectories and utilize a geometric definition of efficiency, η = g φ s (34)Another important observation to be made here isthat, as noted in [2, 3], for the fully actuated swimmers,the gaits that optimized for forward velocity and effi-ciency were the same within each system. There are tworeasons for this:1. At any given power level P , maximum forwardspeed is attained by executing the efficiency-maximizing gait at a pacing that has an even powerusage instead of surging and sagging.2. The shape of the efficiency-maximizing gait is in-variant across power levels because the displace-ment resulting from, and the cost incurred by, exe-cuting one cycle of any gait at optimal pacing doesnot depend on the power level P .The maximum power available, P , dictates the time pe-riod T of the speed-maximizing gait. The attainablemaximum speed increases linearly with the amount ofpower available. For efficiency-maximizing gait, the ef-ficiency is invariant across our choice of time period, T ,because the displacement resulting from and the cost in-curred by executing one cycle of any gait at optimal pac-ing, does not depend on the power level P . However, themaximum power available P sets a lower bound on thepossible range for T .Moving from an active to a passive swimmer also af-fects the gradient calculation process that forms the back-bone of the framework presented in [8]. In a passiveswimmer, the movement of the passive joint is coupledto the movement of the active joint, and therefore wecannot parametrize the motion of the passive joint inde-pendently of the motion of the active joint. We discussthe implications of a passive joint on gait parametrizationand gradient calculation in § IV B. A. Efficiency for swimmers with passive elasticjoints In fully actuated swimmers, the transformation of ef-ficiency from the inverse of cost of transport as definedin (33) to a more geometric definition in (34) was possi-ble because we knew that the optimal pacing keeps therate of power dissipation constant [7]. In the case ofswimmers with passive joints, the response of the passivejoint is dictated by the dynamics of the active joint, asillustrated by the Bode plots shown in Fig. 7. As a re-sult, there is a unique pacing associated with every gaitthe passive swimmer can execute. Changing the pacingof the actuated joint changes the response of the passivejoint as shown in § III and hence the shape of the gait.Thus to find the most efficient gait for the Purcell swim-mer with a passive elastic joint, we have to directly usethe definition of efficiency in (33).While we could choose a constant power pacing forall gaits for fully actuated swimmers, in the case of thePurcell and T-link swimmers with a passive elastic joint,every gait that respects the passive dynamics of the elas-tic joint comes with an inherent pacing. Thus for a givenspring stiffness of the passive joint, the gait that maxi-mizes forward velocity comes with a power requirementassociated with it. Even if we are capable of giving thesystem more power, there is no way for the system toutilize that power to go faster. Hence for the swimmingsystems considered in this paper, there are two mean-ingful measures for comparing different gaits that leadto different definitions of gait optimality: Gaits can becompared by1. Comparing the average speeds they produce ( η = g φ T )2. Comparing their energetic efficiency ( η in (33)) B. Gait parametrization for passive swimmers We use a truncated Fourier series to parametrize thegaits. This choice of parametrization lets us accuratelyapproximate a large family of smooth periodic gaits. Theframework introduced in [8] uses a gradient descent algo-rithm to identify gaits that maximize efficiency as definedin (34). During the gradient calculation process outlinedin appendix B, it is useful to think of the gait as beingparametrized by a series of waypoints. We can generatethese waypoints from the Fourier parametrization. Weuse the gradients calculated at each of these waypointsto calculate gradients with respect to the Fourier seriesparametrization.In the case of swimmers with a passive joint, we letthe actuated joint trajectory α ( t ) be given by a fourthorder Fourier series, α ( t ) = a + (cid:88) i =1 a i cos (cid:16) πiT t (cid:17) + b i sin (cid:16) πiT t (cid:17) (35)Using (25), i.e. the transfer function relating the move-ment of the active and passive joints, we obtain the re- sponse of passive joint to α ( t ) as α ( t ) = L − ( H ( s ) L ( α ( t )) . (36)Using explicit evaluation of the transfer functionfrom (27) and (29), we can write the steady state re-sponse of the passive joint as α ( t ) = (cid:88) i =1 c i cos (cid:16) πiT t (cid:17) + d i sin (cid:16) πiT t (cid:17) , (37)where c i and d i are functions of a i and b i and T .Using this low-order Fourier series parameterization ofthe gait, we can generate the direct transcription way-points, calculate the gradient of the objective function ateach waypoint (details of this gradient calculation processfor speed-maximizing and efficiency-maximizing gaits arepresented in § V and § VI respectively), then project thesegradients onto the Fourier basis, to obtain gradients withrespect to the Fourier series parameters.If the system were fully actuated, we could move theFourier series parameters along these calculated gradientdirections to obtain the optimal gait. In the case of pas-sive swimmers, the Fourier coefficients of the unactuatedshape direction ( c i , d i ) are functions of the Fourier shapecoefficients of actuated shape direction ( a i , b i ). There-fore, to find the correct gradient directions for the Fouriercoefficients of the actuated shape, we have to accountfor the change in the unactuated shape direction that achange in the actuated shape direction would produce.For an objective function f that depends on a i , b i , c i and d i , we can calculate the total derivatives of f withrespect to a i and b i as dfda i = ∂f∂a i + ∂c i ∂a i ∂f∂c i + ∂d i ∂a i ∂f∂d i (38) dfdb i = ∂f∂b i + ∂c i ∂b i ∂f∂c i + ∂d i ∂b i ∂f∂d i , (39)where ∂c i ∂a i , ∂c i ∂b i , ∂d i ∂a i and ∂d i ∂b i are directly taken from thetransfer function coefficients (32) ∂c i ∂a i = A (40) ∂c i ∂b i = A (41) ∂d i ∂a i = − A (42) ∂d i ∂b i = A . (43)We use these total derivatives to calculate the correctgradient directions for the Fourier coefficients of the ac-tuated shape variables, which account for the fact thatin passive swimmers, a change in the shape of input tothe actuated shape variable affects the response of thepassive elastic joint.0 Speed-maximization Efficiency-maximizationFourier series parametriazationof input to actuated joint.(a ,a ,...,a ,b ,...,b ) Time period of input (T)Displacement per cycle( ) Energy per cycle (E)Parameters describing input to actuated jointsOptimality criteria FIG. 8: A flowchart describing how the parameters describing ourinput to the actuated joint affect the two optimality criteria inthe case of drag-dominated swimmers with a passive joint. V. SPEED-MAXIMIZING GAITS A major goal of this paper is to geometrically identifygaits for the Purcell and T-link swimmers with an elas-tic joint that will give it the maximum forward velocity.Therefore the objective function we set out to maximizeis g φ T , where g φ is the displacement in the body x direc-tion produced on executing the gait φ , and T is the timeperiod required to execute the gait.As explained in § IV.B, we parametrize the actuatedshape variable using a low-order Fourier series and ob-tain the Fourier series parametrization of the resultingmotion of the passive joint using the dynamics of thepassive elastic joint presented in § III. With these Fourier-series parameters, p f , we can obtain a sequence of way-points p i (in our implementation we use 100 waypoints),equally spaced in time, that describe the location of thediscretization points in the shape space. As illustrated inthis section, we can then calculate the gradient of speedwith respect to each of these waypoints, i.e., calculate theeffect moving the waypoints would have on the forwardspeed attained by the swimmer on executing the gait. Wethen project these gradients onto the Fourier series ba-sis as explained in § IV.B to obtain the speed-maximizinggaits using gradient descent.We start from the basic variational principle that func-tions reach their extrema when their derivatives go tozero. Given a gait parametrization defined by waypointparameters p , executed in time T , the maximum-velocitycycle must satisfy the condition that the gradient of for-ward velocity with respect to the parameters p and T iszero, i.e. ∂∂p (cid:16) g φ T (cid:17) = 1 T ∂g φ ∂p − g φ T ∂T∂p (44)= 1 T ∂g φ ∂p = 0 , (45) and ∂∂T (cid:16) g φ T (cid:17) = 1 T ∂g φ ∂T − g φ T = 0 , (46)where ∂T∂p term in (44) can be taken as zero since T and p are two independent variables describing the gait, i.e.the time taken to complete a gait T does not depend onthe shape and pacing of the input to the actuated jointdescribed by p .Once we have the gradient of speed with respect to adirect transcription parametrization p obtained from aFourier series parametrization p f , we follow the processoutlined in § IV B, specifically (38) and (39), to obtainthe gradient of speed with respect to the Fourier seriesparametrization, ∂∂p f (cid:16) g φ T (cid:17) . A graphical depiction of howelements of the Fourier series parametrization affect thegait optimization process is shown in Fig. 8.For suitable seed values p f and T , the maximum-velocity gaits can thus be obtained by finding the equi-librium of the dynamical system˙ p f = ∂∂p f (cid:16) g φ T (cid:17) (47)˙ T = ∂∂T (cid:16) g φ T (cid:17) = 1 T ∂g φ ∂T − g φ T . (48)In the geometric framework we introduced in [2], weshowed that the process of finding efficient gaits for thefully actuated Purcell swimmer is akin to the dynamicsof a soap-bubble in which internal pressure and surfacetension combine to determine the shape, size and surfaceconcentration of the soap bubble. This is not the case inswimmers with a passive elastic joint. Instead of two in-dependent processes, the optimization process for findingthe fastest gait is more unified where (B3) is the equationthat helps obtain the shape of the optimal input to theactuated joint and (48) is the equation that helps obtainthe optimal pacing of the input. A. Shape gradient of the optimal input to theactuated joint. As discussed in [8], the gradient that affects the shapeof the input to the actuated joint, dg φ dp , pushes the gait to-wards maximum displacement cycles. From (3), and thefact that variations in p affect φ but not the underlying D A structure, we can reduce this term to ∂g φ ∂p ≈ ∇ p ¨ φ a ( − D A ) . (49)A powerful geometric principle (the general form of the Leibniz integral rule [40]) tells us that the gradient ofan integral with respect to variations of its boundary isequal to the gradient of the boundary with respect to thesevariations, multiplied by the integrand evaluated along theboundary, and allows us to rewrite (49) in terms of the1constraint curvature value and how changes in parame-ter values move the gait’s trajectory through the shapespace.Formally, the multiplication of the gradient of theboundary and the integrand evaluated along the bound-ary is the interior product [41] of the boundary gradientwith the integrand, ∇ p ¨ φ a D (– A ) = ‰ φ ( ∇ p φ ) ⌟ D (– A ) . (50)In systems with just two shape variables, the interiorproduct reduces to a simple multiplication between theoutward component of ∇ p φ and the scalar magnitude ofthe Lie bracket, ∇ p ¨ φ a ( − D A ) = ˛ φ ( ∇ p ⊥ φ )( − D A ) . (51)This gradient calculation is illustrated in Fig. 15. Thegradient of the enclosed area with respect to variations inthe position of p i , i.e. ∇ p i φ a in the e (cid:107) and e ⊥ directions,is the change in triangle’s area as p i moves. Becausethe triangle’s area is always one half base times height(regardless of its pitch or the ratio of its sidelengths),this gradient evaluates to ∇ p i φ a = (cid:2) e (cid:107) e ⊥ (cid:3) (cid:20) (cid:96)/ (cid:21) . (52)Note that this term matches the right-hand side of (51),with only normal motions of the boundary affecting theenclosed area. B. Frequency gradient of the optimal input to theactuated joint In the case of the fully actuated Purcell swimmer, theshape of the gait, and therefore the displacement pro-duced by executing the gait, are independent of the timetaken to execute the gait. This is not true in the case ofthe Purcell swimmer with a passive elastic joint.In this subsection, we examine the gradient that guidesthe optimizer towards the optimal frequency of input tothe actuated joint. When the time period required toexecute the gait is changed, the shape of the gait changesdue to the coupling between the frequency of input tothe actuated joint and the response of the passive jointas described in § III. Changing the time period T thuschanges not only the frequency of the gait cycle but alsothe displacement produced per cycle.We use the chain rule to calculate this gradient, ∂∂T (cid:16) g φ T (cid:17) = 1 T ∂g φ ∂T − g φ T (53)= 1 T (cid:32) ∂g φ ∂α ∂α ∂T + ∂g φ ∂α ∂α ∂T (cid:33) − g φ T (54) Because α is the actuated shape variable and the shapeof the input actuation is independent of the frequency ofactuation, ∂α ∂T reduces to zero. Therefore, the gradientof speed with respect to T reduces to ∂∂T (cid:16) g φ T (cid:17) = 1 T (cid:32) ∂g φ ∂α ∂α ∂T (cid:33) − g φ T . (55)The first term of the right hand side of (55) capturesthe contribution to the velocity of the gait caused by thechange in the shape of the gait resulting from a changein T . The second term accounts for the fact that, evenwithout a change in the shape of the gait, an increase inthe time required to execute the gait would result in adecrease in the velocity of the gait. C. Passive Purcell and T-Link swimmers We implemented the optimizer described in the § V inMatlab by providing (47) and (48) as the gradient for the fmincon optimizer using the sqp algorithm. The shapeof the gait obtained is illustrated in Fig. 9(a). Fig. 9(b)shows the power input to the actuated joint over thecycle. Fig. 9(c) shows a comparison of the speeds achiev-able by the Passive and fully actuated Purcell swimmersat different power levels. Fig. 10 shows the same resultsfor a passive T-link swimmer.The transfer function relating the response of the pas-sive joint to oscillations of the input joint is given by (26), H ( s ) = sC eq s + K eq . Therefore a change in the value ofthe spring stiffness does not affect the fundamental shapeand nature of the response shown by the Bode plot inFig. 7, but it does shift the entire bode plot to the left orright along the frequency axis. Thus an increase in springstiffness shifts the Bode plot to the right, which resultsin the shape of the speed maximizing gait remaining thesame, but the time period required to complete the gaitdecreases, leading to faster speeds. VI. ENERGY-EFFICIENT GAITS In this section, we describe the gradient calculationsinvolved in identifying the gait that maximizes the effi-ciency of the swimmers. The objective function we setout to maximize is η = g φ E , (56)where g φ is the displacement produced on executing thegait φ and E is the total energy expended by the robotexecuting the gait, i.e. E = P avg T (57)2 Units== length of the swimmer= stiffness of the = drag force per unit lengthpassive jointper velocity (c) Maximum speeds achieved by the passive and active Purcell swimmer (b) Shape and power of input to the actuated joint of the passive Purcell swimmer (a) Speed-maximizing gait and pacing FIG. 9: Gaits that maximize speed along x-direction for the Purcell swimmers. (a) The optimal gait for the passive Purcell swimmer(red) and the optimal gait for the fully actuated Purcell swimmer (black). The thickness of the line shows the magnitude of powerrequired at different points of the gait. (b) The power required to execute the optimal gait for the passive Purcell swimmer and input tothe actuated joint over one gait cycle. (c) A comparison table of the speeds achieved by the passive and fully actuated Purcell swimmersat different power levels. (c) Maximum speeds achieved by the passive and active T-link swimmer (b) Shape and power of input to the actuated joint of the passive T-link swimmer (a) Speed-maximizing gait and pacing Units== length of the swimmer= stiffness of the = drag force per unit lengthpassive jointper velocity FIG. 10: Gaits that maximize speed along x-direction for the T-link swimmers. (a) The optimal gait for the passive T-link swimmer(red) and the optimal gait for the fully actuated T-link swimmer (black). The thickness of the line corresponds to the magnitude ofpower required at different points of the gait. (b) The power required to execute the optimal gait for the passive T-link swimmer andinput to the actuated joint over one gait cycle. (c) A comparison table of the speeds achieved by the passive and fully actuated T-linkswimmers at different power levels. As explained in § IV.B, we parametrize the actuatedshape variable using a low-order Fourier series, p f , andobtain the Fourier series parametrization of the result-ing motion of the passive joint using the dynamics of thepassive elastic joint presented in § III. With these Fourier-series parameters, we can obtain a sequence of waypoints p i equally spaced in time, that explicitly define the loca-tion of the discretization points in the shape space. Thesewaypoints form a direct-transcription parametrization ofthe gait p , obtained from the Fourier series parametriza-tion p f . In this section, we present the calculation of thegradient of efficiency at each of these points with respectto p . We then project these gradients onto the Fourierseries basis as explained in § IV B to obtain the efficiencymaximizing gaits using gradient descent.The maximum-efficiency cycle must satisfy the condi-tion that the gradient of efficiency with respect to the parameters, p and T is zero, i.e. ∂∂p (cid:16) g φ E (cid:1) = 1 E ∂g φ ∂p − g φ E ∂E∂p = 0 (58)and ∂∂T (cid:16) g φ E (cid:1) = 1 E ∂g φ ∂T − g φ E ∂E∂T = 0 (59)Once we have the gradient of efficiency with respect toa direct transcription parametrization p obtained from aFourier series parametrization p f , we follow the processoutlined in § IV B to obtain the gradient of efficiency withrespect to the Fourier series parametrization, ∂∂p f (cid:16) g φ E (cid:17) .Thus for suitable seed values of p f and T , the maxi-mum efficiency gaits can be obtained by finding the equi-librium of the dynamical system,3˙ p f = ∂∂p f (cid:16) g φ E (cid:17) (60)˙ T = 1 E ∂g φ ∂T − g φ E ∂E∂T (61)Thus the process of finding the most efficient gait is theresult of a unified process where (60) is the equation thathelps find the shape of the input to the actuated jointand (61) is the equation that helps find the optimal pac-ing of the input. Note that as in § V, the two equationsdo not operate independently, and the gradient of shapedepends on the time period T , and the gradient of thetime period depends on the shape of the gait. A. Shape gradient of the optimal input to theactuated joint The shape of the optimal input to the actuated joint isaffected by two gradients, ∂g φ ∂p and ∂E∂p (60). The details ofhow ∂g φ ∂p is calculated are explained in § V.A. Whereas ∂g φ ∂p pushes the gait towards maximum displacement cycles, E is a measure of the cost required to execute the gaitand ∂E∂p pushes the gait towards low cost shapes. At themost efficient gait, these two opposing gradients canceleach other out, and we get an equilibrium for the gaitoptimization process.Over a gait cycle, no energy is stored in the spring.Hence we can calculate the energy expended, P , whileexecuting a gait by integrating the power flow throughthe actuated joint ( α ). E = ˆ T ˙ α ( t ) T τ dt (62)= ˆ T ˙ α ( t ) T M ( t ) ˙ α ( t ) dt (63)where M ( t ) is the first row of the power metric M ( t ).The gradient of cost with respect to the shape of the gait, ∂E∂p , is calculated by ∂E∂p = ∂∂p ˆ T ˙ α ( t ) T M ( t ) ˙ α ( t ) dt (64)= ˆ T (cid:16) ∂ ˙ α ∂p M ˙ α + (65)˙ α T M ∂ ˙ α∂p + ˙ α T ∂ M ∂p ˙ α (cid:17) dt (66) B. Frequency gradient of the optimal input to theactuated joint The equation that governs the optimization process forfinding the time period of the most efficient gait is de- scribed by (61). The term ∂g φ ∂T is calculated as describedin § V.B The second gradient in the right hand side of (61)is calculated as dEdT = ∂E∂α ∂α ∂T + ∂E∂α ∂α ∂T + ∂E∂T . (67)Because α is the actuated shape variable and the shapeof the input actuation is independent of the frequency ofactuation, ∂α ∂T reduces to zero. Therefore the gradient ofenergy with respect to period reduces to dEdT = ∂E∂α ∂α ∂T + ∂E∂T . (68)The first term accounts for the fact that a change infrequency would change the response of the passive joint α resulting in a change in the shape of the gait, andhence a change in the power dissipated while executingthe gait. The second term accounts for the fact that evenif the shape of the gait remains unchanged, a change inthe frequency of input to the actuated joint will changethe time required to execute the gait and hence wouldchange the power dissipated while executing the gait. C. Passive Purcell and T-Link swimmers We implemented the optimizer described in the § VIin Matlab by providing the gradients of efficiency withrespect to shape and time period, calculated using (60)and (61) respectively, to the fmincon optimizer using thesqp algorithm. The shape of the gait obtained is illus-trated in Fig. 11(a) for a Purcell swimmer. Fig. 11(b)shows the power input to the actuated joint over thecycle. Fig. 9(c) shows a comparison of the efficien-cies achievable by the Passive and fully actuated Purcellswimmers when the forward speed for all the systems isfixed to be equal to the forward speed achieved by thepassive swimmer when executing its maximum efficiencycycle. Fig. 12 shows the same results for a passive T-link swimmer. Note that for the Purcell swimmer, theoptimizer stops because reducing the frequency furtheror making the gait smaller does not provide any mean-ingful increase in efficiency. This observation is line withthe results from [9], where the efficiency was found toasymptotically approach a maximum value as frequencyof gait oscillations approached zero. The maximum effi-ciency gait for the passive T-link swimmer is much largercompared to the maximum efficiency gait of the passivePurcell swimmer because a larger gait helps exploit thepresence of two peaks in the constraint curvature func-tion of the T-link swimmer.As discussed in § V.C, change in spring stiffness doesnot affect the shape and nature of the response of thepassive joint, but shifts the bode plot shown in Fig. 7 tothe left or right along the frequency axis. An increasein spring stiffness shifts the bode plot to the right, whichresults in the shape of the efficiency maximizing gait cycleremaining the same, but the time taken to execute the4 Fully Actuated Purcell Swimmer Passive gait same speed 1.33Ac � ve gait same speed 1.99 Passive Purcell Swimmer Units== length of the swimmer= stiffness of the passive joint (c) Maximum efficiencies achieved by the passive and active Purcell swimmer (b) Shape and power of input to the actuated joint of the passive Purcell swimmer (a) Efficiency-maximizing gait and pacing FIG. 11: Gaits that maximize efficiency along x-direction for the Purcell swimmers. (a) The optimal gait for the passive Purcellswimmer (red) and the optimal gait for the fully actuated Purcell swimmer (black). The thickness of the line shows the magnitude ofpower required at different points of the gait. (b) The power required to execute the optimal gait for the passive Purcell swimmer andinput to the actuated joint over one gait cycle. (c) A comparison table of the efficiencies achieved by the passive and fully actuatedPurcell swimmers moving forward at the same speed as the passive swimmer. Note that for the passive Purcell swimmer, the optimizerstops because reducing the frequency further or making the gait smaller does not provide any meaningful increase in efficiency. Thisobservation is line with the results in [9]. (c) Maximum efficiencies achieved by the passive and active T-link swimmer Units== length of the swimmer= stiffness of the passive joint (b) Shape and power of input to the actuated joint of the passive T-link swimmer (a) Efficiency-maximizing gait and pacing FIG. 12: Gaits that maximize efficiency along x-direction for the T-link swimmers. (a) The optimal gait for the passive T-link swimmer(red) and the optimal gait for the fully actuated T-link swimmer (black). The thickness of the line corresponds to the magnitude ofpower required at different points of the gait. The maximum efficiency gait of the passive T-link swimmer is much larger compared tothe maximum efficiency gait of the passive Purcell swimmer because a larger gait helps exploit the presence of two peaks in theconstraint curvature function of the T-link swimmer. (b) The power required to execute the optimal gait for the passive T-link swimmerand input to the actuated joint over one gait cycle. (c) A comparison table of the efficiencies achieved by the passive and fully actuatedT-link swimmers moving forward at the same speed as the passive swimmer. gait will decrease. The energy required to execute a gaitis inversely proportional to the time taken to execute thegait. Therefore, more energy is required to execute a gaitfaster. Thus, increasing spring stiffness will result in anoverall decrease in the efficiency of swimming. VII. COMPARISON WITH PREVIOUS WORK Passive swimmers have been previously studied, e.g.,in [9–13]. Of these works, [10, 11] and [9] are most rele- vant to this paper, as they discuss the motion of swim-mers with a harmonically-driven active joint and a pas-sive joint. The T-link swimmer used in this paper was in-troduced in [10] as a simplified model to study the swim-ming dynamics of Schistosoma mansoni . The analysisin [9] is particularly relevant, as it applies perturbationtheory to investigate the motion of the Purcell swimmer.The approximation used in the perturbation analysisin [9] is equivalent to assuming that the displacement pro-duced by executing a gait is equal to the area enclosed bythe gait in the shape space multiplied by the constraint5 (c) Purcell Swimmer(d) T-Link Swimmer (g) Purcell Swimmer(h) T-Link SwimmerSpeed vs Frequency Efficiency vs Frequency Efficiency vs Stroke Amplitude(e) Purcell Swimmer(f) T-Link SwimmerSpeed vs Stroke Amplitude(a) Purcell Swimmer(b) T-Link Swimmer Optimal InputSinusoidal Input FIG. 13: In all the subfigures, the solid red lines and the red circles show the speeds and efficiencies predicted by numerical simulationand integral of CCF respectively. The solid black lines show the speeds and efficiencies predicted by the constant-CCF assumption usedin [9] in the link-attached coordinates. The first column of figures illustrate the speed of the (a) Purcell swimmer as a function ofactuation frequency when the input to the controlled joint is a sinusoidal oscillation of unit amplitude and when the input to thecontrolled joint has the optimal shape and (b) T-link swimmer as a function of actuation frequency when the input to the controlledjoint has the optimal shape as obtained in § V. In (a), the grey line shows the speed predicted by the constant-CCF assumption for asinusoidal input. The small amplitude perturbation analysis in [9] predicts speed to be a monotonically increasing function of frequencyfor all inputs. The speed is a monotonically increasing function of frequency for the sinusoidal input. However, it is not monotonicallyincreasing for the optimal input contrary to the prediction from the small amplitude perturbation analysis in [9]. The second column offigures illustrate the efficiency of the (c) Purcell swimmer and (d) T-link swimmer as a function of actuation frequency when the input tothe controlled joint has the optimal shape as obtained in § VI. Figures (e)-(h) illustrate the speed and efficiency of the Purcell swimmerand T-link swimmer as a function of the magnitude of optimally shaped input actuation obtained from § VI. curvature value at the center of the gait, g φ ≈ ¨ φ a D ( − A ) | = φ a · D ( − A ) | . (69)This approximation has been used in several works fromthe geometric mechanics community to identify usefulshape oscillations that resulted in useful net displace-ments e.g., [42, 43].Using the approximation in (69), the authors of [9] con-cluded that for harmonic inputs, the speed of the swim-mer monotonically increases with frequency and asymp-totically approaches a maximum value as the actuationfrequency approaches ∞ . They similarly concluded thatthe efficiency of the swimmer asymptotically approachesa maximum value as the actuation frequency approacheszero.A drawback of this approximation is its the accuracyfalls steeply with increasing gait size, as a larger gaitwould mean larger variations in the value of CCF insidethe region bounded by the gait. Our analysis improveson this approximation in three ways: first, our use of abody-averaged frame instead of the link-attached frame in [9] “flattens” some of the nonlinearity in the system,expanding the domain of gait amplitude for which theperturbation analysis gives accurate results. Second, ouruse of the full integral of constraint curvature over thearea enclosed by the gait, g φ ≈ ¨ φ a D ( − A ) , (70)absorbs much of the remaining nonlinearity. Third, weuse (70) only to calculate gradients, and numerically eval-uate the value of g φ at each step of the gradient descentprocess described in § V and § VI to avoid compoundingerrors from any residual nonlinearityFig. 13(a) and Fig. 13(b) show how the effect on swim-ming speed from changing the frequency of the inputstroke for Purcell and T-link swimmers respectively. Thesolid red lines and the red circles show the speeds pre-dicted by numerical simulation and area integral of cur-vature respectively. The solid black and grey lines showthe speeds predicted by the constant-CCF assumptionused in [9] in the link-attached coordinates for optimallyshaped and sinusoidal inputs respectively. The speedspredicted by the constant-CCF assumption are higher6than the actual speeds obtained by numerical simula-tions, and our integral of CCF is a good approximationof the ground-truth simulation.We can also see from Fig. 13(a) that the velocity ob-tained by the Purcell swimmer does not monotonicallyincrease with frequency for all inputs to the actuatedjoint. This results in the speed-maximizing gaits foundin § V having an optimal frequency associated with them,rather than exhibiting a monotonic increase in speed withfrequency.In the case of efficiency-maximizing gaits, for the op-timal gait shape obtained in § VI, the efficiency doesasymptotically approach a maximum value as shown inFig. 13(c) and Fig. 13(d) as the frequency approacheszero. However, the value is different from the maximumefficiency predicted by applying the small perturbationanalysis from [9] to the T-link swimmer, showing thatsmall perturbation analysis does not completely charac-terize optimal performance. In the case of the Purcellswimmer, the efficiency-maximizing gait is small enoughfor the perturbation analysis to yield accurate results.Figs. 13(e)-(h) illustrate how the constant-CCF as-sumption can introduce errors in identifying optimal ac-tuation shape. Figs. 13(e)-(f) illustrate the effect onthe swimming speed from changing the size of the inputstroke (the reference input is the optimal input obtainedin § V). Figs. 13(g)-(h) illustrate the effect on the swim-ming efficiency from changing the size of the input stroke(the reference input is the optimal input obtained in § VI).The solid red lines and the red circles show the speedsand efficiencies predicted by numerical simulation andintegral of CCF respectively. The solid black lines showthe speeds and efficiencies predicted by the constant-CCFassumption used in [9] in the link-attached coordinates.We can see that in the case of the Purcell swimmer, theconstant-CCF assumption used in [9] incorrectly predictsa monotonic increase in speed with an increase in the am-plitude of the input to the actuated joint. In the case ofthe T-link swimmer, the constant-CCF assumption usedin [9] incorrectly predicts an increase in efficiency as weshrink the optimal gait. The efficiency would go down ifwe shrink the optimal gait, because the CCF value for T-link swimmer is higher at the edges than at the center ofthe shape space. When we shrink the gait, it loses theseregions of high value, leading to a decrease in efficiency,which is not captured by the constant-CCF assumption.The T-link swimmer was first introduced in [10]. Theanalysis in this paper agrees with the most relevant re-sults from [10], which are that when the actuated jointis driven by a simple harmonic input:1. There exists a linear relationship between thespeed-maximizing value of spring stiffness and fre-quency of actuation.2. The average swimming speed increases monoton-ically as the amplitude of actuation is increasedfrom π radians to π radians. In § V.C, we noted that an increase or decrease in springstiffness shifts the Bode plot of the response of the pas-sive joint to the right or to the left without changingthe shape of the Bode plot resulting in frequency of thespeed-maximizing input being linearly related to the stiff-ness of the passive joint.In Fig. 10(a), we can see that the CCF value for T-link swimmer is higher at the edges than at the centerof the shape space. Thus an increase in the amplitudeof actuation would enclose more of the high value regionleading to an increase in speed. VIII. CONCLUSIONS In this paper, we have identified the geometric struc-ture of optimal gaits for viscous swimmers with passiveelastic joints by combining the constraint-curvature anal-ysis in [8] with frequency-response models for the steadystate motion of driven oscillators. We use this struc-ture to identify both speed-maximizing and efficiency-maximizing gaits. The optimal gaits for passive swim-mers maximize the CCF integral relative to perimeterand pacing costs, subject to amplitude and phase con-straints of a first order system.As discussed in § IV, for the fully actuated swimmers,the maximum forward speed achievable is only restrictedby the maximum power we are able to supply the joints,but for the swimmers with the passive elastic joint, evenwith more powerful actuators, there is a theoretical max-imum forward speed the system can achieve dictated bythe stiffness of the passive joint.The important factor that makes the performance ofthe fully actuated swimmers superior to that of the swim-mers with the passive elastic joint in terms of energy ef-ficiency is the fact that not only can the fully actuatedswimmers execute a much larger set of gaits, they canexecute any gait the passive swimmer can execute at apacing that is just as good or better than the pacingdictated by the dynamics of the passive joint.This raises the question of what benefits, e.g., simplic-ity of construction, does having a passive elastic mem-ber give to biological organisms that locomote in a lowReynolds number fluid? Most biological organisms havetails that resemble an elastic filament. The propulsiveand flexive dynamics of such filaments have been wellstudied [1, 44–46]. Artificial microscopic swimmers withelastic filaments have been proposed based on this bodyof work [47, 48]. An interesting line of future work wouldinvolve investigating the tradeoff between elastic elementinefficiencies and structural complexity of being fully ac-tuated.This work is the first step towards expanding the appli-cability of the geometric framework presented in [2, 3, 8]to systems where underactuated shape parameters playa role in the dynamics of the system. In the case of thepassive Purcell swimmer, assuming the torque required toaffect a desired shape change did not depend on the cur-7rent shape of the swimmer did not introduce significanterrors in the predicition of the limit cycle correspondingto inputs to the actuated joint as shown in Fig. 4. Thismight not always be the case in all the swimmers we con-sider. A future line of research would be to improve ourfrequency domain analysis by using non-linear perturba-tion theory to obtain more accurate predictions of limitcycles.Another line of future work would involve studying theshape of optimal gaits in swimmers with more than onepassive joint (e.g. a four-link swimmer with two passivejoints). The relationship between design choices (e.g. ra-tio of link lengths, ratio of joint stiffness) and the shapeof the optimal gaits would also be an interesting questionto answer in systems with more than one passive joint. Appendix A: Comparison with Lighthill Efficiency Lighthill’s efficiency, defined as the (reciprocal) ratiobetween the average power consumed by a given strokeand the power that would be required to drag the swim-mer at the same average velocity as that produced by thestroke, η lh = P ref P avg , (A1)is a commonly-used measure of swimming performance.Because the drag force acting on the swimmer is pro-portional to the velocity of its motion, it is readily shownthat for a given swimmer, Lighthill’s efficiency is propor-tional to v P avg [9]. Ignoring the proportionality factor, wecan thus take the Lighthill efficiency as η lh = v P avg , (A2)where v avg is the average velocity of the swimmer.In [8], we use a geometric measure of efficiency, η ,defined in (34) as the ratio of displacement produced percycle to the pathlength of the cycle in the shapespace(weighted by the shapespace metric M ), η = g φ s (A3)to identify optimal gaits for fully actuated swimmers.In this appendix, we demonstrate that for fully actu-ated swimmers the gait that maximizes η , maximizesLighthill’s efficiency and vice versa.Because s is equal to the time integral of the squareroot of instantaneous power expended, we can rewrite η as η = g φ ´ T (cid:112) P ( t ) dt . (A4) Surface tension from distance metricOutward pressure from Lie bracket FIG. 14: Our algorithm in [8] maximizes gait efficiency in fullyactuated swimmers by finding cycles in the space of body shapesthat enclose the most curvature of the system dynamics(measured via the Lie bracket) while minimizing theircost-to-execute (measured as the metric-weighted lengths of theirperimeters). This process is analogous to the process by which airpressure and surface tension combine to produce the shape andsize of a soap bubble. Top: The forward progress of a locomotingsystem as it executes a gait cycle. From [7], we know that the optimal pacing for any gaitutilizes a constant power pacing and hence for all t , P ( t ) = P avg (A5)Substituting (A5) into the expression for η provides η = g φ (cid:112) P avg T (A6)= v avg (cid:112) P avg (A7)= √ η lh (A8)Because the square root is a monotonic function, a gaitthat maximizes our definition of efficiency with respectto path and is executed at constant power pacing alsomaximizes the Lighthill efficiency optimized over pathand pacing, and vice versa.Note that this conclusion does not hold for the under-actuated systems we consider in this paper, for whichpath and pacing cannot be controlled individually. Appendix B: Soap-bubble Gait Optimization In this appendix, we present a brief overview of theframework introduced in [8] to identify optimal gaitsfor fully-actuated drag-dominated swimmers. In § Vand § VI, we build on this framework to identify efficiency-maximizing and speed-maximizing gaits for passive8 Gait described byparameters p. Moving in the direction producesno change in area Change in area produced due tomovement along direction Movement in the direction does not change position of along the axis Movement in the direction balances the side lengths Change in metric and segment length with movement of point FIG. 15: Changes in area caused by moving in the two coordinate directions in the local frame. Moving in the tangential direction e (cid:107) produces no change in area, as the area of the triangle given by half the product of base length and height remains the same. swimmers with one passive and one active joint. Theframework uses a gradient descent algorithm to identifygaits that maximize efficiency as defined in (34).Given a gait parametrization p, maximum-efficiencycycles satisfy the condition that the gradient of the effi-ciency ratio is zero, ∇ p g φ s = 1 s ∇ p g φ − g φ s ∇ p s = . (B1)For suitable seed values p , solutions to (B1) can there-fore be reached by finding the equilibrium of the dynam-ical system ˙ p = ∇ p g φ s . (B2)The stable equilibria of the right-hand equation in (B2)are gaits in the same “image families” as the system’soptimally-efficient gaits (i.e., they follow the same curveas the optimal gait, but not necessarily at the samepacing). To construct the optimal gait, we can eitheroptimize via (B2) and then choose a constant-metric-speed parameterization, such that the pacing penalty σ from (17) goes to zero, or directly include ∇ p σ in ouroptimizer[49].Combining the gradient of the pacing term with thegradient of the image-optimizer places the maximum-efficiency gait as the equilibrium of˙ p = ∇ p g φ − g φ s ∇ p s + ∇ p σ (B3)(from which we have factored out a coefficient of s from (B1)).As illustrated in Fig. 14, this differential equation is di-rectly analogous to the equations governing the shape ofa soap bubble: ∇ p g φ takes the Lie bracket as an “internal pressure” seeking to expand the gait cycle to fully encir-cle a sign-definite region, ∇ p s is the “surface tension”that constrains the growth of the bubble, and ∇ p σ is the“concentration gradient” that spreads the soap over thebubble’s surface.For the fully actuated swimmers, we parametrize thegait as a sequence of waypoints p i such that the gait pa-rameters p i explicitly define the location of the discretiza-tion points. As illustrated in Fig. 15, each waypoint p i forms a triangle with its neighboring points and we candefine a local tangent direction e (cid:107) as p i +1 − p i − = (cid:96) e (cid:107) (B4)and a local normal direction e ⊥ orthogonal to e (cid:107) .We selected this direct-transcription parameterizationbecause it facilitated visualizing the workings of our op-timizer (and thus the dynamics governing any other op-timization applied to this problem). Additionally, it al-lowed us to illustrate simultaneous optimization of thegait path and its pacing. We could also parametrize thegait using a Fourier series or Legendre polynomials. Inthis case, the pacing optimization should be done afterthe image of the optimal gait has been found becausefinding an optimal pacing can no longer be formulated asa process orthogonal to the gradient descent process forfinding the image of the optimal gait. ACKNOWLEDGEMENT This work was supported by the National ScienceFoundation, under CMMI grants 1462555 and 1653220.The authors would like to thank Shai Revzen andBrain Bittner for several productive discussions that con-9tributed to this work, and Saad Bhamla for bringing thework in [10] to our attention. [1] C. Brennen and H. Winet, Fluid mechanics of propulsionby cilia and flagella, Annual Review of Fluid Mechanics , 339 (1977).[2] S. Ramasamy and R. L. Hatton, Soap-bubble optimiza-tion of gaits, in Decision and Control (CDC), 2016 IEEE55th Conference on (IEEE, 2016) pp. 1056–1062.[3] S. Ramasamy and R. L. Hatton, Geometric gait opti-mization beyond two dimensions, in American ControlConference (ACC), 2017 (IEEE, 2017) pp. 642–648.[4] R. L. Hatton and H. Choset, Geometric swimming atlow and high Reynolds numbers, IEEE Transactions onRobotics , 615 (2013).[5] R. L. Hatton and H. Choset, Nonconservativity and non-commutativity in locomotion, European Physical JournalSpecial Topics: Dynamics of Animal Systems , 3141(2015).[6] R. L. Hatton, T. Dear, and H. Choset, Kinematic car-tography and the efficiency of viscous swimming, IEEETransactions on Robotics (2017).[7] L. Becker, S. A. Koehler, and H. A. Stone, On self-propulsion of micro-machines at low Reynolds number:Purcell’s three-link swimmer, Journal of Fluid Mechan-ics , 15 (2003).[8] S. Ramasamy and R. L. Hatton, The geometry of opti-mal gaits for drag-dominated kinematic systems, IEEETransactions on Robotics , 1014 (2019).[9] E. Passov and Y. Or, Dynamics of purcells three-linkmicroswimmer with a passive elastic tail, The EuropeanPhysical Journal E , 1 (2012).[10] D. Krishnamurthy, G. Katsikis, A. Bhargava, andM. Prakash, Schistosoma mansoni cercariae swim effi-ciently by exploiting an elastohydrodynamic coupling,Nature Physics , 266 (2017).[11] A. Montino and A. DeSimone, Three-sphere low-reynolds-number swimmer with a passive elastic arm,The European Physical Journal E , 42 (2015).[12] I. Jo, Y. Huang, W. Zimmermann, and E. Kanso, Passiveswimming in viscous oscillatory flows, Physical Review E , 063116 (2016).[13] L. J. Burton, R. L. Hatton, H. Choset, and A. E. Hosoi,Two-link swimming using buoyant orientation, Physicsof Fluids , 091703 (2010).[14] E. M. Purcell, Life at low Reynolds numbers, AmericanJournal of Physics , 3 (1977).[15] In the parlance of geometric mechanics, this assigns Q the structure of a (trivial, principal) fiber bundle , with G the fiber space and R the base space .[16] This kinematic condition has been demonstrated fora wide variety of physical systems, including thosewhose behavior is dictated by conservation of momen-tum [19, 22], non-holonomic constraints such as passivewheels [20, 22, 42, 50], and fluid interactions at the ex-tremes of low [4, 23] and high [4, 21, 51] Reynolds num-bers. [17] R. L. Hatton and H. Choset, Geometric motion plan-ning: The local connection, Stokes’ theorem, and the im-portance of coordinate choice, International Journal ofRobotics Research , 988 (2011).[18] S. D. K. Kelly and R. M. Murray, Geometric phases androbotic locomotion, J. Robotic Systems , 417 (1995).[19] G. C. Walsh and S. Sastry, On reorienting linked rigidbodies using internal motions, Robotics and Automation,IEEE Transactions on , 139 (1995).[20] J. P. Ostrowski and J. Burdick, The mechanics and con-trol of undulatory locomotion, International Journal ofRobotics Research , 683 (1998).[21] J. B. Melli, C. W. Rowley, and D. S. Rufat, Motion plan-ning for an articulated body in a perfect planar fluid,SIAM Journal of Applied Dynamical Systems , 650(2006).[22] E. A. Shammas, H. Choset, and A. A. Rizzi, Geometricmotion planning analysis for two classes of underactuatedmechanical systems, Int. J. of Robotics Research , 1043(2007).[23] J. E. Avron and O. Raz, A geometric theory of swim-ming: Purcell’s swimmer and its symmetrized cousin,New Journal of Physics (2008).[24] This approximation (a generalized form of Stokes’ the-orem) is a truncation of the Baker-Campbell-Hausdorfseries for path-ordered exponentiation on a noncommu-tative group, and closely related to the Magnus expan-sion [52, 53]. For a discussion of the accuracy of this ap-proximation and its derivation, see [5].[25] R. L. Hatton and H. Choset, Optimizing coordinatechoice for locomoting systems, in Proceedings of the IEEEInternational Conference on Robotics and Automation (Anchorage, AK USA, 2010) pp. 4493–4498.[26] This model is most widely associated with swimmersat low Reynolds numbers (e.g., [34]),but can also beregarded as an informative general model for systemsthat experience more lateral drag than longitudinal drag(e.g. [54]). Our choice of resistive force here also doesnot preclude the use of more dejointed physical models(e.g., [55]) to construct the local connection A .[27] Note that a more complete fluid model for low Reynoldsnumber swimming including interbody flow interactionswould change numbers in the dynamics but not the over-all structure.[28] A constraint that the allowable velocities are orthogonalto a set of locally-linear constraints, i.e., that they are inthe nullspace of a constraint matrix ω .[29] This metric corresponds to approximating the viscousdrag acting on the swimmer via a local resistive forcemodel. If we used a more accurate fluid model, the dragon a low-Reynolds number swimmer would have the sameform, except that C would also depend on r and (cid:96) .[30] J. E. Avron, O. Gat, and O. Kenneth, Optimal swimmingat low reynolds numbers, Physical Review Letters ,186001 (2004). [31] R. L. Hatton and H. Choset, Kinematic cartography andthe efficiency of viscous swimming, IEEE Transactionson Robotics , 523 (2017).[32] A. E. Bryson and Y.-C. Ho, Applied Optimal Control,Optimization, Estimation, and Control (John Wiley &Sons, New York-London-Sydney-Toronto, 1975) pp. 402–402.[33] J. P. Ostrowski, J. P. Desai, and V. Kumar, Optimal gaitselection for nonholonomic locomotion systems, Interna-tional Journal of Robotics Research , 225 (2000).[34] D. Tam and A. E. Hosoi, Optimal stroke patterns forPurcell’s three-link swimmer, Phys. Review Letters ,068105 (2007).[35] A. DeSimone, L. Heltai, F. Alouges, and A. Lefebvre-Lepot, Computing optimal strokes for low Reynolds num-ber swimmers, Natural Locomotion in Fluids and on Sur-faces , 177 (2012).[36] L. Giraldi, P. Martinon, and M. Zoppello, Controllabilityand optimal strokes for n-link microswimmer, in Decisionand Control (CDC), 2013 IEEE 52nd Annual Conferenceon (2013) pp. 3870–3875.[37] P. Bettiol, B. Bonnard, L. Giraldi, P. Martinon, andJ. Rouot, The Purcell Three-link swimmer: some geo-metric and numerical aspects related to periodic optimalcontrols (2015), working paper or preprint.[38] A. Shapere and F. Wilczek, Geometry of self-propulsionat low Reynolds number, Journal of Fluid Mechanics , 557 (1989).[39] J. Lighthill, Flagellar hydrodynamics, SIAM Review ,161 (1976).[40] H. Flanders, Differentiation under the integral sign, TheAmerican Mathematical Monthly , 615 (1973).[41] Not the inner product, see [40] for more details.[42] R. M. Murray and S. S. Sastry, Nonholonomic motionplanning: Steering using sinusoids, IEEE Transactionson Automatic Control , 700 (1993).[43] K. A. Morgansen, B. I. Triplett, and D. J. Klein, Geomet-ric methods for modeling and control of free-swimmingfin-actuated underwater vehicles, IEEE Transactions onRobotics , 1184 (2007).[44] C. H. Wiggins and R. E. Goldstein, Flexive and propul-sive dynamics of elastica at low reynolds number, Physi-cal Review Letters , 3879 (1998).[45] R. E. Goldstein and S. A. Langer, Nonlinear dynamics ofstiff polymers, Physical review letters , 1094 (1995). [46] T. S. Yu, E. Lauga, and A. Hosoi, Experimental investi-gations of elastic tail propulsion at low reynolds number,Physics of Fluids , 091701 (2006).[47] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A.Stone, and J. Bibette, Microscopic artificial swimmers,Nature , 862 (2005).[48] J. Edd, S. Payen, B. Rubinsky, M. L. Stoller, andM. Sitti, Biomimetic propulsion for a swimming surgi-cal micro-robot, in Proceedings 2003 IEEE/RSJ Inter-national Conference on Intelligent Robots and Systems(IROS 2003)(Cat. No. 03CH37453) , Vol. 3 (IEEE, 2003)pp. 2583–2588.[49] Including ∇ p σ in the optimizer works best for parame-terizations in which ∇ p σ is orthogonal to ∇ p g φ s , such aswaypoint based direct transcriptions. For other parame-terizations, e.g., Fourier series, the gradients may not beorthogonal and a two-step procedure of optimizing theimage then the pacing will produce better results. Forwaypoint-based parameterizations, the ∇ p σ term has asecondary benefit of helping to stabilize the optimizer bymaintaining an even spacing of points, and thereby pre-venting the formation singularities in the curve).[50] A. M. Bloch et al. , Nonholonomic Mechanics and Control (Springer, 2003).[51] E. Kanso, Swimming due to transverse shape deforma-tions, Journal of Fluid Mechanics , 127 (2009).[52] J. E. Radford and J. W. Burdick, Local motion plan-ning for nonholonomic control systems evolving on prin-cipal bundles, in Proceedings of the International Sympo-sium on Mathematical Theory of Networks and Systems (Padova, Italy, 1998).[53] W. Magnus, On the exponential solution of differentialequations for a linear operator, Communications on Pureand Applied Mathematics VII , 649 (1954).[54] R. L. Hatton, H. Choset, Y. Ding, and D. I. Goldman,Geometric visualization of self-propulsion in a complexmedium, Physical Review Letters , 078101 (2013).[55] N. Giuliani, L. Heltai, and A. DeSimone, Predicting andoptimizing microswimmer performance from the hydro-dynamics of its components: The relevance of interac-tions, Soft robotics5