Computational Design of Skinned Quad-Robots
Xudong Feng, Jiafeng Liu, Huamin Wang, Yin Yang, Hujun Bao, Bernd Bickel, Weiwei Xu
AAUGUST 2019 1
Computational Design of Skinned Quad-Robots
Xudong Feng, Jiafeng Liu, Huamin Wang,
Member, IEEE
Yin Yang,
Member, IEEE
Hujun Bao,Bernd Bickel, and Weiwei Xu,
Member, IEEE
Abstract —We present a computational design system that assists users to model, optimize, and fabricate quad-robots with soft skins.Oursystem addresses the challenging task of predicting their physical behavior by fully integrating the multibody dynamics of the mechanicalskeleton and the elastic behavior of the soft skin. The developed motion control strategy uses an alternating optimization scheme to avoidexpensive full space time-optimization, interleaving space-time optimization for the skeleton and frame-by-frame optimization for the fulldynamics. The output are motor torques to drive the robot to achieve a user prescribed motion trajectory.We also provide a collection ofconvenient engineering tools and empirical manufacturing guidance to support the fabrication of the designed quad-robot. We validatethe feasibility of designs generated with our system through physics simulations and with a physically-fabricated prototype.
Index Terms —Computational Fabrication, Motion Design, 3D-Printing, Physics-based Simulation. (cid:70)
NTRODUCTION T H e design and construction of robots that can execute com-pelling motions is a challenging task. It requires careful geo-metric planning of robotic mechanisms and professional knowledgeof the kinematic and dynamic behavior of the robot. Embeddingsuch knowledge into procedures of computational robot design [1],[2] in conjunction with rapid prototyping techniques, such as3D printing technology, bears tremendous potential to acceleratethe construction of personalized robots. For instance, Megaro etal. [3] used a kinematic optimization algorithm for the design ofmultilegged robots consisting of rigid links. However, real-worldcreatures are not merely rigid skeleton rigs. The muscle and fleshsurrounding the skeleton provides diverse morphologies, enrichedexpressivity. For instance, people wearing a prosthetic limb oftenprefer a highly realistic rubbery artificial arm over a more functionalmechanical one [4]. Moreover, skins and muscles might be essentialfor facilitating challenging tasks, such as reproducing compliantgrasping of a hand or swimming motions of a fish [5], [6].Different to the simulation of the robots with rigid links only,the influence of the soft body on the control torques at actuatorsand on the contact forces at the ground should be simulated andfully taken care of to judge the plausibility of a designed soft robot.This imposes a significant computational challenge for the motiondesign and fabrication of skinned robots. In this paper, our goalis to address this problem by integrating both dynamic simulationand kinematic optimization into the motion design of quad-robotsystems, with soft skins attached as their organic embodiments.Its kernel is an optimization problem that integrates both user-provided kinematic preference and physical constraints of the robotto obtain a dynamically feasible motion plan. The primary physicalconstraint is the dynamic viability of the skeleton trajectory when • Xudong Feng, Jiafeng Liu, Hujun Bao and Weiwei Xu are with State KeyLab of CAD& CG, Department of Computer Science, Zhejiang university,Zhejiang 310058, China. Weiwei Xu is the corresponding author.E-mail: { bao, xww } @cad.zju.edu.cn • Huamin Wang is with Department of Computer Science and Engineering,Ohio State University. E-mail:[email protected]. • Yin Yang are with the Department of Electrical and Computer Engineering,University of New Mexico. E-mail: [email protected]. • Bernd Bickel are with Institute of Science and Technology Austria, Vienna,Austria. E-mail: [email protected]. a soft skin is attached. The dynamics of the robot is formulated astwo-way coupled subsystems of the rigid multibody system and thedeformable skin. In addition, our system also incorporates motorconstraints and a stability constraint. The motor constraints ensurethat the calculated joint torques are within the physical limits ofthe installed servo motors, whereas the stability constraint requiresthe center of projection (COP) of the robot structure to fall insidethe supporting polygon. As a result, given the surface mesh anddesired kinematic trajectories of the robot, our pipeline generates aphysically-valid motion plan that can be realized under the givenhardware constraints. The robot itself can be fabricated using rapidprototyping technology.While space-time optimization is widely used in long-horizonmotion design problem, it is computationally prohibitive if directlyapplied to our case of two-way coupled system, mostly due tothe large number of DOFs of the soft skin mesh coupled with theskeleton and various physcial constraints. To this end, we proposean alternating optimization algorithm of two optimization stepsto mitigate the computational cost: (1) Space-time optimizationwith repect to the DOFs of the rigid skeleton while assuming thedeformation of the soft skin remains the same as the previousiteration. (2) Frame-by-frame optimization as in [7] to furtheroptimize the motion plan obtained in step 1 according to all thephysical constraints at each frame. To this end, we non-triviallyextend the spring-based control force formulation in [7] to handlethe full simulation of the two-way coupled rigid skeleton and softskin. Our solver can efficiently handle the Lagrange multipliersintroduced by the coupling constraints and the collisions betweenthe skeleton and skin. With a proper initialization, these two stepsare alternatively executed until the convergence. In the space-timeoptimization step, the influence of the soft skin to the rigid skeletonis treated as known quantity and simplified to be the coupling forcesand the influence of the center of mass (COM) at each link in theskeleton due to the skin deformation. This setting makes the space-time optimization computationally efficient through decouplingthe skin mesh DOFs, which is critical to achieve global effect inthe motion design. The skin deformation is updated after eachframe-by-frame optimization.To ease fabrication, we provide a convenient workflow withtailored engineering tools and empirical fabrication guidance a r X i v : . [ c s . G R ] J u l UGUST 2019 2
Input mesh Mechanical skeleton & skin mesh Kinematic trajectories Motion design resultFabrication result Walking robot
Fig. 1. We propose a computational fabrication system for designing and fabricating skinned quad-robots. Given an input mesh representing theshape of a quad-robot such as this beetle-like robot, we design its mechanical skeleton with motors to drive its locomotion. The motion plan isgenerated by an optimization algorithm with kinematic trajectories as the user input. The trajectories consist of the foot swing trajectories (red curves),center of mass (COM) trajectory (blue curve) and foot contact plan (yellow bars) input by the user. Our optimization fully takes account of all thephysical and dynamical constraints. By fabricating this robot design, it can be verified that our algorithm is able to generate plausible and physicallyfeasible motion plans for quad-robots, and the simulated results well match the physical experimental results. Note that we only render the transparentinput surface without thickness to demonstrate the surface-structure coupling geometry. embedded in a standard CAD software, empowering regular usersto design quad-robots. The modular design scheme allows the userto quickly start from a design template of the mechanical skeletonin
SolidWorks and adapt it to the body shape. The rigid skeletonis fabricated via 3D printing, and the skin is separately fabricatedusing injection molding by pieces. We tested the optimizationalgorithm on skinned quad-robots with varying body-to-leg ratiosand different mechanical skeletons. Both physical and numericalexperiments show that the proposed algorithm is an effective meansof obtaining physically valid motions of the skinned quad-robots.
ELATED W ORK
Computational fabrication aims at designing and creating physi-cal artifacts with the help of computational methods. A large classof methods addresses inverse design problems by incorporating fab-rication limitations in geometric design algorithms via constrainedoptimization or the integration of fast simulation techniques [8],[9]. This line of research enables the design of objects with a widerange of controllable physical and mechanical properties, such asappearance [10], [11], [12], [13], deformation [14], [15], [16], [17],articulation [18], [19], [20], and mechanical motion [21], [22], [23],[24], [25]. Some existing contributions also investigated how toinstantiate virtual characters as 3D-printable physical entities likemechanical robots [3], [21]. Yet, these methods merely focus onrobots consisting of rigid links and basic balancing constraintsand/or velocity limits.Bickel et al. [26] proposed a process for designing syntheticskin and actuation parameters for animatronic characters that mimicfacial expressions of a given subject. Skouras et al. [15] optimizedthe internal material distribution so that the resulting characterexhibits the desired deformation behavior. Focusing on actuation,Bern et al. [27] computed the layout of winch-tendon networksto animate plush toys, and Ma et al. [28] optimized the chamberstructure and material distribution for designing soft pneumaticobjects.Our work shares some of these goals but takes a significantlydifferent approach. Instead of relying on a rigid or quasi-static underlying simulation of the skin deformation in our case, theinfluence of the soft skin on the embedded moving mechanicalskeleton makes us face a dynamic two-way coupled multibody-elastic problem. It is much more challenging to accurately simulateand optimize due to the drastically increased system complexity,nonlinearity and discontinuity, i.e., due to the complimentaryconstraints enforced at the contact vertices.
Physics-based character motion generation has vast applicationsin both graphics and robotics. Algorithmic approaches include lever-aging space-time optimization with necessary physical constraintsand developing controllers to drive forward simulations.The seminal work by Witkin and Kass [29] generated motiontrajectories by optimizing physical constraints and animator con-trols at key frames, a well-known space-time constraints frameworkfor animation. With proper motion data, the space-time optimizationproduces realistic articulated motions for bipedal or multileggedcharacters through different physical properties [30], [31], [32],[33], [34], [35], [36], [37], [38]. It can be used to transform motioncapture data into physically plausible motions [39].The locomotion controller aims to compute joint torques orcontrol forces to drive the locomotion behaviors of articulatedfigures. The joint torques are usually calculated via the proportionaland derivative (PD) controller such that the rigid skeleton of acharacter follows designated joint angle trajectories [40], [41], [42].Balance control strategies, such as the swing foot placement or zeromoment point, are essential to generating stable locomotions. [40],[41], [43], [44], [45], [46], [47], [48]. Continuous adaptation of thetarget joint trajectory for balancing a walking human was developedin [49]. Controllers that produce highly dynamic skills for humananimation were suggested in [50], [51], [52]. The joint torques canalso be computed via optimal control to approximate the motioncapture data or motion data from kinematic simulators [53], [54].Our work is inspired by studies on how to drive the soft skindeformation with the underyling rigid skeletons or pseudo muscleforce [7], [55], [56]. Two-way coupling of rigid bodies and elasticbodies was considered in [57]. Fast simulation and control ofsoft robots of various configurations and actuations has also been
UGUST 2019 3 studied using finite element method and the reduced formulationof compliance matrix [58], [59], [60], [61], [62].
Elastic body simulation focuses on the formulation of an elasticdeformation energy and the proper handling of contact constraintsto simulate realistic deformations of soft bodies [63], [64], [65],[66]. A comprehensive survey of physics-based elastic deformationmodels can be found in [67].Space-time optimization techniques can also be applied to con-trol the motion of elastic bodies that are represented by volumetricmeshes. To reduce the number of variables used to control thevertex positions in the optimization, model reduction techniquesare frequently used [68], [69], [70]. Barbiˇc et al. [68] imposed theequation of motion constraint in elastic body deformation, usingthe discrete adjoint method to compute the gradients of controlforces. Pan et al. [71] integrated the contact forces as additionalvariables to handle environment interactions and solved the space-time objective with alternating optimization, but did not handle thetwo-coupling problem we want to solve.
VERVIEW
Given an input surface mesh of a quad-robot, we first designmechanical skeleton and skin mesh (Section 6.1) for the robot, andthen use the proposed two-step alternating algorithm to optimizefor a physically plausible motion plan (Section 5). The overallsystem flowchart is illustrated in Fig. 1.During the alternating optimization, the first space-time op-timization step outputs joint angle trajectories for the designaccording to the user-specified end-effector and COM trajectories(Section 5.1). This step is made possible by only consideringthe approximated skin deformation. The second frame-by-frameoptimization step improves the physical plausibility of the jointangle trajectories with full simulation and various physical con-straints(Section 5.2), such as physical torque limits for the selectedmotors in the design. Two-way coupled multibody-elastic dynamics(Section 4) are adopted in this step for the full simulation.Finally, the designed robot is fabricated by fast prototypingmethods for a physical validation (Section 6.2), and stepper motorsare mounted to drive the skeleton and the attached soft skin torealize the motion plan. WO - WAY C OUPLED M ULTIBODY -E LASTIC D Y - NAMICS
A core ingredient of our system is the dynamic simulation ofthe robot using a self-actuated rigid skeleton with a soft skinattached. We are inspired by existing coupled simulation systemsin graphics [72], [73], [74] and exploit the Lagrange multipliersmethod to enforce the two-way coupling between the skeleton andthe soft skin, which can be naturally integrated into the subsequentlocomotion optimization. The Lagrange multipliers are used toguarantee that the skeleton and skin are attached to each other atprescribed locations, and we solve for all the unknown DOFs fromboth subsystems simultaneously.
We use Lagrangian mechanics for both the rigid skeleton andsoft skin and obtain a symmetric formulation for these twosubsystems. To ease the formulation of motion contraints, we usethe generalized coordinates q to parameterize the entire skeleton, where q = { c x , c y , c z , q , ..., q m } . The vector q is composed of theCartesian coordinates of the center of mass (COM) at the rootlink { c x , c y , c z } and the three Euler angles at each joint { q , ..., q m } .We opt to model the soft skin using the neo-Hookean materialmodel [75] because it has been demonstrated to be well suited forrobotic skins made out of silicone [76] and can handle large localdeformations induced by joint rotation observed in our examples.If necessary, however, our approach should be easily extensibleto more sophisticated material models such as Mooney-Rivlinand Ogden. The DOFs representing the vertex displacements ofthe tetrahedral skin mesh are denoted by the vector u . Thoughnumerical discretization, the equations of the forward simulationof rigid skeleton and soft skin at each time step can be writteninto a linear system A x = b , where A is the system matrix and b is a vector of the sum of constant terms and external forces. Thedetailed derivation of two linear systems, A r , b r for rigid skeletonand A d , b d for soft skin, are elaborated in Appendix A.The two subsystems are coupled by attaching the soft skin tothe underlying skeleton at prescribed glue vertices . Mathematically,this straightforward treatment leads to a set of nonlinear positionconstraints: C ( q , u ) = R ( q ) r + t − x c = . (1)The matrix R converts the positions on the rigid links where theskin is attached, r , from local to world coordinates. This operationcan be easily expressed as a rotation chain from the links all the wayback to the root. Meanwhile, t concatenates the root translations ofthe rigid links, and x c denotes the positions of the glue vertices onthe skin mesh x c = S c u , where S c is a selection matrix.Using the Lagrange multipliers method, we obtain the coupledmultibody-elastic system: A r ∇ q C (cid:62) d ∇ u C (cid:62) ∇ q C ∇ u C ∆ q ∆ u λλλ = b r b d . (2)The coupling constraint is linearized via ∇ C . In each time step,we solve for the changes of the system DOFs in Eq. 30. Thus, wehave q i = q i − + ∆ q and u i = u i − + ∆ u , where the superscript [ · ] i indicates the frame index. There are two types of collisions/contacts we need to take careof in our simulation. The first is the collision between the robot’sfeet and the ground surface, which provides necessary supportand friction forces to the robot. As will be detailed in Sec. 5.2,we resolve them using linear complementary constraints (LCP) toguarantee the physical feasibility of the locomotion.The second is the self-collision of the soft skin; rotatingjoints tend to compress the inward skin and induce self-collisions.Moreover, skin-skeleton collisions could also occur when theskeleton is being articulated. Because such collisions typicallytake place under low relative velocities, we handle them using theexplicit penalty force method [77], [78]. OTION P LAN O PTIMIZATION
As the system input, the trajectories of the COM, end effectors,and the footfall pattern are provided by the user. Our system
1. As discussed in Sec. 6, skin-skeleton collisions can occur between the skinand the bracing unit, not the mechanical skeleton.
UGUST 2019 4 generates a dynamically feasible motion plan for the robot thatresembles as much as possible the one prescribed by the user.A motion plan, defined as P = { q i , i = ... N ; ∆ t } , includes theskeleton configuration q i of the i -th frame for i = N , and thetime step size ∆ t .Previous works, e.g. [Megaro et al 15], solved the motiondesign problem by simultaneously finding the optimal skeletonconfiguration for all the frames in P . However, this problembecomes much more challenging in the case of a multibody-elastic system, due to the large number of DoFs and associatedphysical constraints. Thus, the global approach with full DOFs isinfeasible in our case. In this section, we elaborate the details ofour two-step alternating motion plan optimization algorithm whichuses approximated skin deformation to significantly improve theefficiency of the global space-time optimization. Specifically, theinfluence of the skin deformation on the skeleton is approximatedas the coupling forces at glue vertices and the influenced COMpositions of each link. The convergence of the algorithm is testedconsidering the difference between the space-time optimized oneand the optimized one after the frame-by-frame optimization step.Fig. 2 illustrates the algorithm flowchart. Similar to the formulation in [36], for each i -th frame, we considerthe set of generalized coordinates of the rigid skeleton, q i , togetherwith the contact forces F ic , j and torques τττ ic , j that are exerted on the j -th end effector in contact with the ground. The influence of theskin deformation to the skeleton at i -th frame is simplified to bethe coupling forces λλλ i at glue vertices and the skin deformation u i .They are obtained from the previous frame-by-frame optimizationand not optimized in this step. The displacements in u i arerepresented into the local coordinate frames of links at each framein order to compute the influence of the skin deformation to theCOM of the robot.Given these quantities, the optimization objective is defined asthe weighted sum of six terms that balance the user-specified endeffector and COM trajectories and the smoothness of the optimizedmotion: E A = min ∑ i (cid:0) α t E i τττ + α s E iS + α c E iCOM + α e E iEE + α f E iF + α o E iO (cid:1) . (3)The first term E i τττ is standard in space-time optimization to minimizethe torques τττ i exerted at the joints: E i τττ = m (cid:13)(cid:13) τττ i ( q i , λλλ i , F ic , τττ ic ) (cid:13)(cid:13) . The torques τττ i are computed using inverse dynamics, and thecoupling forces λλλ i are integrated as the external forces exerted bythe elastic skin at the coupling points.The second term E iS encourages the smoothness of the opti-mized motion, which is: E iS = (cid:13)(cid:13) q i + − q i + q i − (cid:13)(cid:13) . The two terms, E iEE and E iCOM , enforce the optimized motionto follow the user-specified end-effector and COM trajectoriesrespectively: E iEE = (cid:13)(cid:13) φ iEE ( q i , u i ) − e i (cid:13)(cid:13) , E iCOM = (cid:13)(cid:13) φ iCOM ( q i , u i ) − g i (cid:13)(cid:13) , Space - time Optimization (Skeleton DOFs)
Initialization
Frame - by - frame Optimization (Full DOFs) q q ’ |q- q’ | End Skin mesh deformation + coupling force No Fig. 2. Alternating algorithm Flowchart. where the functions φ iEE ( q i , u i ) and φ iCOM ( q i , u i ) define how tocompute the end effector and COM positions, given the generalizedcoordinates of the skeleton and the deformation of the skin mesh.For each end effector, we select one vertex closest to the endeffector of the rigid skeleton and represent this vertex into the localcoordinate system of the end effector to compute φ EE (see Fig. 3).The term E iF penalizes the deviation from the motion { ˜ q i , i = , .., N } generated in the previous frame-by-frame optimization orthe initialized motion plan at the first iteration: E iF = (cid:13)(cid:13) q i − ˜ q i (cid:13)(cid:13) . After the frame-by-frame optimization step, the weight of the E F term to follow the motion plan of the previous iteration in thespace-time optimization is increased by 10%, which means thealgorithm leans toward physical plausibility.The last term E O enforces that the end-effectors in the contactare flat: E iO = (cid:13)(cid:13) ψ iEE ( q i ) − ˆ n (cid:13)(cid:13) , where ψ iEE ( q i ) is a function to compute the orientation of the endeffector and ˆ n is set to be ( , , ) , the normal of the support plane.The formulation of this term is motivated by the fact our skinnedrobot is soft and hence a contact area appears whenever feet contactthe ground. We try to maximize the contact area at the moment ofcontact because it is important to achieve a stable motion. Similarly,when the foot is about to hover, we want to clear as many contactvertices as possible. Therefore, when an end effector is close tothese important moments, we add E O to the objective function,which try to make the end effector face the upright direction of theground.We also impose hard kinematic and dynamic constraints on theoptimized variables to obtain a stable motion plan. The kinematicconstraints include the contact constraint, the center of pressure(COP) constraint and the optional periodic constraint, while thedynamic constraints include the momentum and the friction forceconstraints as in [36]. Contact constraint : In the optimization, we need to enforce thefootfall pattern that is specified by the user to encode when the footshould leave or touch the ground. This constraint can be writteninto: c ij φ iEE , j ( q i , u i ) y = , ∀ i , j , c i − j c ij ( φ iEE , j ( q i , u i ) − φ i − EE , j ( q i , u i )) = , ∀ i , j , (4)where c ij is a binary variable set to 1 if the j -th end effector isin contact with the ground at the i -th frame. This variable canbe directly derived from the footfall pattern. The y coordinate ofthe end effector, denoted by φ iEE , j ( q i , u i ) y , should be 0 since theground is set to be y = COP constraint : The COP should be inside the supporting polygonfor a stable motion plan. This constraint can be written into: P · φ iCOP ≤ , ∀ i , (5) UGUST 2019 5 Fig. 3. Foot contacts. where the funtion φ iCOP computes the COPposition using the same method as in [3].The rows in P represent edges of thesupporting polygon. Since the space-timeoptimization requires the foot to be flaton the supporting plane, the supportingpolygon is formed by the convex hull ofthe vertices that represent the sole meshesof the end-effectors in contact. However,the contact forces and torques at a sole aresimplified to be exerted on a single pointof the end effector to ease the optimization. Fig. 3 illustrates thecontact points (red balls) selected as the positions of the end-effectfor the beetle-like robot. Periodic constraint : When the user desires a periodic motion, thejoint angles are expected to be the same in the first and the lastframe of the optimized motion plan: J ( q ) = J ( q N ) , (6)where J extracts the joint angles from the generalized coordinates. Momentum constraint : The change of the linear and angularmomentum of the robot should be determined by the externalcontact forces and torques, which can be formulated into thefollowing equations:˙ R i = mg + ∑ j c ij F ic , j , ∀ i , j , ˙ L i = ∑ j c ij (cid:0)(cid:0) p ic , j − φ iCOM ( q i , u i ) (cid:1) × F ic , j + τττ ic , j (cid:1) , ∀ i , j , (7)where R i and L i are, respectively, the total linear and angularmomentum of the robot at i -th frame, and p ic , j gives the contactposition of the j -th end effector. Friction force constraint : The contact force should be inside thefriction cone to satisfy the Coulomb model of friction. Specifically,we have: 1 m ( µ F ic , j ⊥ − (cid:107) F ic , j (cid:107) (cid:107) ) ≥ , ∀ c ij = , m ( (cid:107) F ic , j (cid:107) − (cid:107) τττ ij (cid:107) ν b (cid:107) − τττ ij ⊥ ν t ) ≥ , ∀ c ij = , (8)where F ic , j ⊥ and F ic , j (cid:107) represent the components of the contactforce perpendicular and parallel to the ground respectively. Thefirst equation requires the contact force to be inside a frictioncone. The second equation relates the contact force and the contacttorque, where ν b and ν t are set to be the radius of the circumcircleof the sole mesh. Optimization : With the defined objective function and constraints,the space-time optimization step can be written into:min q i , F ic , j , τττ ic , j , i = ,.. N E A st. : φ e = , φ g ≤ φ e and φ g represent the equality and inequality constraintsrespectively. We solve this sequential quadratic programmingproblem using the Gauss-Newton algorithm. The weights in theobjective function are specified as follows: α t = e − , α s = . , α c = α e = α f = α o = In the frame-by-frame optimization, we drop the simplificationof the previous step and consider the full dynamics formulated in Eq. (30) in a similar way as in [7]. This allows us to furtherimprove the physical plausibility of the motion. In practice, thisrequires the solution of a difficult quadratic programming problemwith complementarity constraints (QPCC) to handle contact andfriction. Different to the soft body only formulation in [7], thecoupling constraints between multi-body skeletion and elastic skinintroduce a large number of Lagrange multiplier variables andlargely increase the complexity of the solver. To speed up thesolution, we follow the condensation technique widely used inphysical simulation [79], [80]. Specifically, we select the drivingtorques at joints and foot contact forces as optimization variablesand lump the DOFs of mesh vertices and rigid skeletons to thesevariables through the condensation of the system matrix in Eq. (30).This choice facilitates the formulation of the physical torque limitconstraints of motors as well. In this section, we describe the detailsof the matrix condensation and the per-frame optimization problem. System matrix condensation : With the coupled multibody-elasticsystem defined in Eq.(30), the nonlinear relation between thesimulated DOFs ∆ u and ∆ q and the vector b r and b d can berevealed by eliminating the unknown Lagrange multipliers λλλ in thematrix condensation (see Appendix.B for the detailed derivation): ∆ u = A − d b d − A − d ∇ u C (cid:62) A − C (cid:0) ∇ q C A − r b r + ∇ u C A − d b d (cid:1) , ∆ q = A − r b r − A − r ∇ q C (cid:62) A − C (cid:0) ∇ q C A − r b r + ∇ u C A − d b d (cid:1) , (10)where A C = ∇ q C A − r ∇ q C (cid:62) + ∇ u C A − d ∇ u C (cid:62) .According to the formulation of b r and b d in Appendix.A, thejoint torques and the contact forces exerted on the soft skin areabsorbed into the vector g r and g d . Since the rest terms in b r and b d are constant, we can use φ ∆ q and φ ∆ u to represent Eq. (10) moreprecisely: ∆ q = φ ∆ q (cid:0) τττ , F ⊥ , F (cid:107) (cid:1) , and ∆ u = φ ∆ u (cid:0) τττ , F ⊥ , F (cid:107) (cid:1) . (11)Here, F ⊥ and F (cid:107) are magnitudes of normal and tangent forcesat all the contact vertices on the skin mesh, and τττ represnets thejoin torques. To handle LCP constraints, the contact force at acontact vertex is modeled as F = ˆ n F ⊥ + D F (cid:107) instead of the 3Dvector representation in Section 5.1, where ˆ n F ⊥ represents theupright supporting force along the contact normal ˆ n = [ , , ] (cid:62) and F ⊥ ∈ R is the force magnitude. D ∈ R × is a matrix and itscolumns are the vectors that span the contact plane. In our system,we use four directions to form the friction cone [81] and F (cid:107) ∈ R is the tangent magnitude. The explicit penalty forces for resolvingthe collision between the skeleton and skin are pre-determinedquantities and merged into the constant terms in b r and b d .Given q i − and u i − , ∆ q i and ∆ u i determine the positions andorientations of COM, COP, and end effectors at the i th frame.Thanks to Eq. (11), these kinematics variables are now functions of τττ , F ⊥ , and F (cid:107) . Therefore, we can derive the functions required inthe computation of kinematic information, namely φ iEE , ψ iEE , φ iCOM and φ iCOP , if external forces τττ i , F i ⊥ , and F i (cid:107) are given. Note thatthe purpose of these functions are explained in the space-timeoptimization (see Sec. 5.1). Optimization : We follow the control strategy used in [7] tooptimize the input motion plan on a frame-by-frame basis. Itis used to make sure that the output joint angle trajectories ofthe space-time optimization step are physically feasible, which isachieved using the two-way coupled multibody-elastic dynamicsas constraints. At each frame, it can be formulated as a quadraticprogramming problem with complementarity constraints.Specifically, we seek for joint torques ( τττ i ) and contact forces( F i ⊥ , F i (cid:107) ) such that the corresponding ∆ q i = φ i ∆ q and ∆ u i = φ i ∆ u UGUST 2019 6 Fig. 4. Our QPCC solver converges quickly in most cases. The left plot isthe converging curve of a frame when the front left leg of the monster-likerobot leaves the ground. The middle plot is the converging curve of aframe when this leg is in the air (i.e. other three feet are on the ground).The right plot is the converging curve of a frame when this leg hits groundagain. satisfy necessary hard constraints and the resulting locomotionmatches the input locomotion as much as possible. Mathematically,it is formulated asmin τττ i , F i ⊥ , F i (cid:107) , λλλ i (cid:107) E G (cid:16) τττ i , F i ⊥ , F i (cid:107) (cid:17) subject to: (cid:107) τττ im (cid:107) < U m , m = , ..., M P · φ iCOP ( τττ i , F i ⊥ , F i (cid:107) ) ≤ ≤ F i ⊥ F i (cid:107) λ i (cid:107) ⊥ ˆ n (cid:62) ∆ u ic ∆ t D (cid:62) ∆ u ic ∆ t + λ (cid:107) µ F i ⊥ − (cid:62) F i (cid:107) ≥ . (12)In Eq. (12), the first hard inequality constraint of (cid:107) τττ im (cid:107) < U m is the motor constraint requiring for all the M motors that thecomputed torque magnitude is within its physical limit U m . Thesecond inequality constraint P · φ COP ≤ requires the positionof the COP to be within the supporting polygon as in Sec. 5.1.The last complementary constraint is enforced at each individualcontact vertex. It characterizes the contact mechanism such thatwhen normal force exists, the relative velocity between the groundand the contact vertex along the contact normal should be zero,etc. Here, ∆ u ic ∈ R is the incremental displacement of a contactvertex. The auxiliary vector λ (cid:107) is related to the tangent velocity ofa sliding contact; µ is the friction coefficient; and is a vector ofones, that is, = [ , , , ] (cid:62) .The objective function E G has four terms: E G = α S E S + α F E F + α O E O + α τττ E τττ + α C E C , (13)where: E S = (cid:13)(cid:13)(cid:13) φ i ∆ q − φ i − ∆ q (cid:13)(cid:13)(cid:13) , E F = (cid:13)(cid:13) q i + φ i ∆ q − ¯ q i + (cid:13)(cid:13) , E O = (cid:13)(cid:13) ψ iEE − ˆ n (cid:13)(cid:13) , E τττ = (cid:107) τττ i − τττ i − | E C = (cid:13)(cid:13) E (cid:0) φ iEE − φ i − EE (cid:1)(cid:13)(cid:13) . (14)The first energy term E S is the smoothness penalty, whichfavors motions with consistent velocities. The term E iF penalizesthe deviation from the motion { ¯ q i , i = , .., N } generated in theprevious space-time optimization. E O is the same soft constraint onthe orientation of an end effector as in Eq. (3), which can maximizethe contact area for a stable motion. E τττ is used to penalize thelarge variation of the control torques at joints between frames.The last term E C imposes a penalty to moving end effectors whoare responsible for supporting feet. In other words, if a foot isin contact with the ground and supporting the body, we use E C to reduce the risk of its possible tangent sliding. Here, E is anelementary matrix that picks positions of supporting end effectors.The weighting constants for each of these penalty terms are asfollows: α S = α F = α O = α τττ = . α C = The key to solving the QPCC problem of Eq. (12) is to have afeasible configuration for all the contact vertices. Our strategy issimilar to that of [7]: We flip complementary constraints whenthe inequality constraint reaches the boundary. Specifically, contactvertices fall into one of the three following categories: • Contact breakage means that the contact vertex will leavethe ground plane in the next frame, and the complementaryconstraints should be lifted. • Sliding indicates that the contact vertex is moving within theground plane. In this situation, the complementary constraint forits normal force F ⊥ becomes: F ⊥ > , ˆ n (cid:62) ∆ u ic ∆ t = , (15)and the complementary constraints for the tangent force F (cid:107) are: F (cid:107) ≥ , D (cid:62) ∆ u ic ∆ t + λ (cid:107) = ; λ (cid:107) ≥ , µ F ⊥ − (cid:62) F (cid:107) = . (16) • Static friction implies that the contact vertex is fixed on thecontact plane. In this case, the complementary constraint for itsnormal force is the same as Eq. (15). The constraints for thetangent force are: F (cid:107) ≥ , D (cid:62) ∆ u ic ∆ t + λ (cid:107) = ; λ (cid:107) = , µ F ⊥ − (cid:62) F (cid:107) ≥ . (17)The inequality constraint of µ F ⊥ − (cid:62) F (cid:107) ≥ τττ i , γγγ c E G (cid:0) τττ i , γγγ c (cid:1) subject to: (cid:107) τττ im (cid:107) < U m , m = , ..., M P · φ iCOP ≤ (cid:2) φ i ∆ u (cid:3) c = , (18)where (cid:2) φ i ∆ u (cid:3) c returns the incremental displacements of all thecontact vertices. This assumption of fixing all the contact verticesis realized via the Lagrange multipliers method, and the resultingmultipliers γγγ c correspond to the constraint forces at these vertices.Now, let γγγ c ∈ R be the constraint force at one of the contactvertices. It can be decomposed along normal and tangent directionsas: γ ⊥ = ˆ n (cid:62) γγγ c , and γγγ (cid:107) = (cid:0) I − ˆ n ˆ n (cid:62) (cid:1) γγγ c . (19)We label all the contact vertices as contact breakage, static friction,or sliding by checking γ ⊥ and γγγ (cid:107) . If γ ⊥ ≤ 0, which indicates acontact breakage, we remove the constraint at the vertex in the nextiteration. If γ ⊥ > 0, we further examine the magnitudes of µγ ⊥ and (cid:107) γγγ (cid:107) (cid:107) . If µγ ⊥ > (cid:107) γγγ (cid:107) (cid:107) , the vertex falls into the static friction category,otherwise the vertex is considered sliding. After all the contact UGUST 2019 7 vertices are labelled, we can convert the complementary constraintsinto a set of equality or inequality constraints, as explained inEqs. (15), (16), and (17), and re-solve the QP optimization. Itis known that QPCC is NP-complete, and few contact verticescould make the optimization procedure computationally intractable.Therefore, we simplify this procedure by grouping vertices on theplanar surface of the foot mesh into five patches similar to [7].In our experiments, we found that such initial vertex groupingoften provide a good start for the QPCC solver. Typical convergingcurves are plotted in Fig. 4, and we stop the optimization after 10iterations. We observe that the condensed QPCC solver is around50x faster than the QPCC without condensation. Given the mechanical skeleton and the skin mesh of a robot, wefirst associate the mesh vertices to the links of the skeleton toobtain its skinning information. Hence, the mesh vertices can bedeformed with the skeleton in the space-time optimization step,while the local coordinates of the vertices should be computedusing their deformed positions and the local frame of the links ateach frame. The initial motion plan are computed through the space-time optimization step without the trajectory following terms E iF .In this step, the skin deformation is assumed to be static and eachlink has additional weights from its associated vertices. Afterwards,the initial skin mesh deformation is simulated by imposing thecoupling constraints in the elastic simulation of the skin, and theinitial coupling force is then obtained according to the deformationof the tetrahedra connected to the coupling points [56]. ESIGN AND F ABRICATION Designing and fabricating a quad-robot is a challenging task.We facilitate the design by using a set of mechanical skeletontemplates, and narrow the gap between professional and regularusers by creating several SolidWorks scripts. This allows even aninexperienced user to tweak high-level semantic parameters. Fig. 6shows three built-in mechanical skeleton templates provided in oursystem for quad-robots. Each template is built of modularized CADparts to ease the fabrication cost. The first one is the design usedin the beetle-like robot, which consists of a torso structure and fourlimb structures. Their exploded views are detailed in the figure aswell. The torso structure has four shoulder joints that connect to itsfour limbs. A microcontroller board sending trigger signals to themotors is mounted inside of the torso. The limb structure includeslinkage parts of an upper leg, a lower leg, and a foot. On each limb,two uniaxial motors are mounted to provide necessary rotationalfreedoms at the knee and the ankle. The other two templates varyin different initial poses and foot link geometries.In the following, we describe the details of the design pipelineand the fabrication procedure respectively. The design starts with a given 3D model that corresponds to theappearance of the robot. Our system extracts an initial skeletalline using the mesh contraction method [82] as shown in Fig. 5.This skeleton is actually an approximation of the medial axisof the model, and it is used as a general guide for the follow-up template embedding and editing. We employed the modulardesign idea so that the user can edit the geometry of a templatemechanical component to obtain a customized mechanical skeleton Torso b m i L Limb structureUpper legLower leg FootKnee AnkleSkin attachmentTorso structureLeft shoulders Right shoulders Fig. 6. Three mechanical structure templates and the exploded views ofthe torso and limb structures of the first template. m m . m m m m m m mm mm L i n k l e n g t h M o t o r b a y Fig. 7. With the assistance of the developed Solidwords scripts, theuser only needs to tweak semantic parameters like the link length, motormount size, etc. to obtain a customized link. The size of all the pilot holesfor screw installation remains unchanged under such edits. for quad-robots of various morphologies. To this end, several SolidWorks scripts are developed to assign semantic parameters Fig. 5. The initialskeletal line of thebeetle-like robot. of a link, such as the link length, motormount size, etc., and the user only needsto tweak these intuitive parameters to ob-tain a personalized design without creatingone from scratch. The size of pilot holeson the link for screw installation remainsunchanged under such edits. An example isgiven in Fig. 7, where the lengths of the linkand the motor bay are increased. Although afew iterations may be necessary during em-bedding, the developed SolidWorks scriptsgreatly accelerate the procedure.Typically, given a new surface model of a quad-robot, weembed limbs first and then adjust the geometry of the torso tomake sure it fits the exterior skin. Specifically, a global scale of themechanical structure template and local rotations of the links arefirst performed so that the template can be inside the input surfacemesh. Then, the user can select the start and end points of a link onthe extracted skeletal line and trigger the designed script to edit thelink geometry to match the specified length and adjust the widthof the link. Finally, the bracing unit of the torso is generated in asimilar way as the skin creation(the details follows shortly), andwe dig out holes to reduce its weight, for instance, the bracing unitfor the torso of the beetle-like robot shown in the first picture inthe second row of Fig. 1. Skin and folding regions creation : The exterior skin of the robotis designed to be 8 mm thick at the foot and 4 mm thick at otherparts by default, and it is created by the mesh hollow operation in UGUST 2019 8 Sweeping contourTeeth width Folding regionCentral axisRotation Fig. 8. We add folding regions to facilitate the stretching deformationof the skin. The template folding region is similar to gear teeth, andit is formed by a sweeping cut operation, that is, the CSG differencebetween the volume surrounded by the original skin surface and thevolume formed by rotating the sweeping contour along a central axis. Materialize Magics . This operation treats the space surroundedby the input 3D surface model as a solid and hollows the interiorspace of the solid to match the specified thickness parametersto create the skin that is amenable to fabrication. When it isbeing bent, the skin can yield rather large resisting forces understretching deformations. Regular commercial motors may notpossess sufficient power to overwhelm the internal stretching. Toresolve this practical issue, we add a few folding regions on theoriginal skin mesh, as shown in Fig. 8. The folding region iscreated by applying a sweeping cut operation in Solidworks overthe original skin surface where the motor is installed. This smalltreatment increases the skin area where substantial bending occursand effectively reduces the resulting stretching force. Note thatthe hollow operation is performed after the creation of the foldingregions. Fortunately, the two software are compatible in mesh fileformat. Fig. 9. Glue vertices. The mechanical structure of the robotis 3D printed with polypropylene-likestereolithography (SLA) resin, which isa widely used material for fabricatingjoints and low-friction moving parts. Theexterior skin of the robot is made of a layerof soft rubbery material and fabricated viainjection molding. To reduce the effort andcost of creating the skin molds, we fabri-cate the skin on a piece-by-piece basis: onelimb has one skin piece, and the torso hastwo pieces as shown in Fig. 1 (Skin pieces).The skin-skeleton attachment is physicallyrealized by another 3D printed bracingunit between the skin and the skeleton.This bracing unit is attached to each linkon the skeleton and serves as a supportingstructure between the rubbery skin and themechanical skeleton (see Fig. 1 & Fig. 9). The purpose of thisdesign is to expect that the friction between the skin and the printedparts can disable the relative motion between skin and skeletonat these parts, which is verified in the physical validation. Wethus select the glue vertices in Fig. 9 according to the position ofskin-skeleton attachment parts so that the coupling constraints canreflect this physical setting. The mass matrix and the inertia tensorof this bracing unit are integrated in our multibody subsystemdynamics. Finally, skin pieces are glued together after all skinpieces are installed using nonreactive PVA adhesive. Motor specification We use the MG995R servo motor to drive themotion of the skinned robot. The motor’s size is 40 . × × 38 mm Robot Joints Skin Glue Weight Opt.Beetle-like 12 14 , 152 140 3 . kg ∼ . s Monster-like 16 16 , 624 301 11 . kg ∼ . s Lizard-like 17 18 , 261 176 10 . kg ∼ . s TABLE 1Physical and simulation statistics of three tested robots. Joints: thenumber of joints on the skeleton. Skin: the number of vertices on the skinmesh. Glue: the number of glue vertices. Weight: the physical weight ofthe robot. Opt.: the average time used for optimizing Eq. (12) of onemotion frame.Fig. 10. We plot the torque values of the optimized motion planning ofthe beetle-like robot. Orange curve: the torque curve at an ankle joint (atthe orange box). Blue curve: the torque curve at a shoulder joint (at theorange box). with the maximum torque of 20 kg · cm under 6 . U m = . · m in Eq. (12)). In total, 12 motors are installed in the beetle-likerobot. All the motors are controlled with an Arduino board, whichsupports up to 32 motors. XPERIMENTAL R ESULTS In this section, we first report the torque limits and folding regionexperiments in the motion design of the beetle-like robot to validateits physical feasibility. The mechanical skeleton of this robot isdesigned based on the first template in Fig. 6 and fabricated using3D printing. A comparison to kinematic optimization only is alsoprovided for this robot. Second, we report the motion design resultsfor two additional quad-robots: a monster-like robot and a lizard-like one. The performance of our optimization algorithm dependson the number of vertices on the skin mesh, the number of gluevertices, and the number of joints of the mechanical skeleton. Theframe interval ∆ t is 0 . 005 second, and we employ the discretecollision detection algorithm to handle self-collisions and foot-ground collisions. Normally a motion cycle has around 500 frames.Our optimization algorithm was implemented on a desktop PCwith an intel i7-7700 CPU and 16 GB memory. The soft skinis made of isotropic rubber material whose Young’s modulus is0 . GPa , and Poisson’s ratio is 0 . 46. Table 1 reports some essentialphysical and simulation statistics of these three examples. Thespace-time optimization step with approximated skin deformationfor a skeleton is around 40 seconds. For the slow-walking motionsas shown in Fig. 1 and Fig. 19, we only need one iteration toconverge, since the frame-by-frame optimization can reproducethe motion from space-time optimization step well with physicalconstraints. For the relatively fast trotting motion in Fig. 18, thealgorithm converges after two iterations. Torque limit : The torque limit in the motion optimization (i.e.Eq. (12)) of the beetle-like robot is set as 1 . 96 N · m to match thephysical torque limit of the used MG995R motor. To verify if thishard constraint is faithfully enforced during the optimization, weexamine torque values that are calculated by our motion design UGUST 2019 9 Without folding With folding Stress distribution Fig. 11. Adding folding regions significantly relieves the stretching stressover the skin. We program the motor at an ankle joint of the beetle-likerobot to rotate ± ◦ within 2 seconds in this experiment. A smooth skincan only be bent around ± ◦ , while the folded skin is able to reach thedesired deformation. The physical experiment (black skins) results areconsistent with simulation results (yellow skins). Frame 390 Frame 936 Ours[Megaro et.al. 2015] Seconds -1.5-1-0.500.5 ee n K -0.5-0.2500.250.5 Knee (Kinematic)Knee (Inveres dynamics)Knee (Our)Root roll (Kinematic)Root roll (Inveres dynamics)Root roll (Our) ll o R Fig. 12. The motion plan generated by the kinematic-only optimization [3]leads to unstable walking sequences. Left: Selected frames. Right: Thejoint angle curves. Kinematic: the kinematic optimization result. Inversedynamics: the simulation result by following the kinematic optimizationresult. Our: our optimization result. The large roll angles of the root link isdue to the unstable pose using kinematic-only optimization. Please seethe accompanying video for the full comparison. system after per-frame optimization. The result is reported inFig. 10, where torque curves at an ankle joint and a shoulderjoint are plotted. It can be seen that the imposed motor constraintsuccessfully bounds the torque magnitude to be within the limit toensure that the designed motion is physically possible. Folding regions Adding folding regions to the robot’s skin isan effective fabrication artisanry to enable the robot assemblyusing off-the-shelf servo motors and lower the fabrication cost.To demonstrate its advantage, we compare the skin deformationunder the joint rotation when the robot is attached to a regularsoft skin and a folded skin. Both skins are fabricated using thesame materials. It can be seen from Fig. 11 that rotating jointsyield large stretching stress over the skin, which could easily gobeyond the physical capacity of many commercial servo motors.In this test, we follow the aforementioned motor specification bysetting the maximal torque to 1 . 96 N · m and test if this power issufficient to generate the necessary skin deformation. We set ourtarget bending angle to ± ◦ , which is a common value in manywalking gaits. The motor is programmed to reach this target in 2seconds. Our simulation shows that the smooth robot skin withoutfolding region prevents the motor from producing sufficient jointrotation and the maximum angle that can be reached is only about ± ◦ . With the folding region, swept by an 8 mm-depth tooth overthe smooth skin, our simulation predicts that the motor is able togenerate the desired rotation. The physical experiment results arequite consistent with our simulation prediction as reported side byside in Fig. 11. VS. kinematic-only optimization : In contrast to robots withonly rigid mechanical skeletons, skinned robots exhibit a muchmore complicated dynamic behavior, which should be fullyincorporated during the motion design. To illustrate the necessity O u r s K i n e m a ti c s on l y Fig. 13. Physical experiments show that kinematic-only optimization isnot a feasible solution for skinned robots – there are noticeable backsteps (highlighted with a red box) in a motion cycle as the driving torques,after damped by the skin deformation, are not strong enough to producenecessary normal contact forces. Please refer to the accompanied videofor a clearer comparison. of incorporating influences of the soft skin, we compare the motionplans generated using our method and the one by Megaro et al. [3].Because the primary focus in [3] is to design robot creatures withonly rigid links, Megaro and colleagues used a kinematic-basedoptimization strategy, which includes the trajectories of COM, COP,end effectors, and the footfall pattern. Based on the resulting motionplan, we compute the corresponding driving torques at joints usinginverse dynamics. Specifically, the driving torques are computedby imposing another set of rotation constraints over the skeletonin Eq. (30) using the Lagrange multipliers method (with necessarycomplimentary constraints and inequality constraints to handlethe ground contact and motor torque limit). The constrained jointrotation corresponds to the one obtained from the kinematic-basedmotion plan, and the multipliers represent required generalizedconstraint force, which are converted to joint torque via J ω k toachieve the target joint rotation. As shown in Fig. 12, the physicalsimulation results suggest that a kinematically valid joint trajectorydoes not guarantee a smooth walking cycle of the skinned robot,even though the constraints of COM/COP are also specified in thekinematic optimization without skin information. The coincidence,for the knee joint, of the joint angle curves of kinematic andinverse dynamics shows that our inverse dynamics computationcan track the kinematic motion plan well, and the roll angle of theroot link experiences a larger variation in the inverse dynamicssimulation. This is also verified in the physical experiment asshown in Fig. 13. The motion plan obtained using only kinematicoptimization leads to a shaky motion. We also observe backwardmotions as highlighted in the figure. Our method, because it fullyconsiders various physics conditions and constraints, yields a muchbetter result. Motion design results : Fig. 14 illustrates the simple foot liftingmotions designed by our system for the beetle-like robot. Thesetwo motions, i.e., single-foot lifting and double-foot lifting, are UGUST 2019 10 (a)(b) Fig. 14. The foot lifting motion for the Beetle-like robot. (a) Single-footlifting. (b) double-foot lifting. The green balls indicate the COP positionsand gray lines the support polygons. Our optimization algorithm canconstrain the COP to be inside the support polygon. Seconds -1.25-1-0.75-0.5-0.25 J o i n t a ng l e -0.05-0.02500.0250.05Knee angles (Kinematic input)Knee angles (With skin weight)Knee angles (With skin deformation)COM Z (Kinematic input)COM Z (With skin weight)COM Z (With skin deformation) C O M Z Fig. 15. The comparison of joint angle curves and COM positions. Withskin weight: the curves are from the initial space-time optimization resultwith only skin weight. With skin deformation: curves from the space-timeoptimization result with skin deformation simulated by frame-by-frameoptimization. COM Z: the z component of the COM, representing theCOM movement from left to right during the motion. also used to show the COP is constrained to be inside thesupport polygon with our optimization algorithm. Please see theaccompanying video for the animation.The embedded skeleton of the monster-like robot shown inFig. 6 is designed using the third template. The weight of this robotis 11 . kg , and its size is 48 . × × 27 cm. The youngmodulus at the tail and belly of the monster is reduced by 85%to demonstrate the dynamics of the skin. Two different input foottrajectories are used to generate the walking motions for the robot.As shown in the leftmost column in Fig. 17, the first trajectory has alonger stride length but lower step, while the second one has shorterstride and higher step. Our system is able to accommodate suchvariations and produces a smooth and physically correct motionplan. The walking speed for these two walking motions are 0 . . 06 meter/second. Frame 25 Frame 148 Fig. 16. The effect of the COP constraints in the trotting motion plan.Red balls: The COP positions computed using the motion plan afterinitialization. Green balls: The COP positions optimized with the skinmesh deformation. The supporting polygons are in cyran. We generate a trotting motion of 0 . 45 meter/second speed forthe monster-like robot to demonstrate that our system can supportfast motion. In this example, the trajectories of joint positionsare labelled using the horse motion pictures photographed byEadweard Muybridge, a famous photographer for his work onmotions. The joint positions are mapped to a horse motion withthe specified speed using the method in [83] and re-targeted tothe skeleton of our monster. This initial kinematic motion (pleasesee the accompanying video for the motion) is then optimizedusing our alternating motion optimization algorithm to turn it intoa physically feasible motion.We notice that the flighting phase of the initial kinematic motionis not consistent with the foot contact plan. To be more specific, thetime of the flighting phase is not enough for the monster to returnback to the ground. Thus, the space-time optimization with physicalconstraints, especially the momentum constraint, is necessary toeliminate such inconsistency, which is critical to the success ofthe alternating optimization. Fig. 16 illustrates the effect of theCOP optimization. The red balls indicate that the COP positionsat two frames in the initialized trotting motion plan are outsideof supporting polygon after the first space-time optimization thatdoes not account for the deformation of the skin mesh. Thus, theframe-by-frame simulation fails to produce a stable trotting motionwith its initially optimized motion plan, as shown in the second rowof Fig. 18. The COP constraint is turned off to produce this failedexample once this constraint can not be satisfied by the solver.After the second iteration, the COPs are moved into supportingpolygons, as indicated by the green balls. A smooth trotting motioncan then be generated as shown in the first row on the right ofFig. 18. The comparisons of the optimized joint angles and the z components of COP are illustrated in Fig. 15. The variation of z components indicates the COP moves from left to right so that it isinside the supporting region.Another example is reported in Fig. 19. The mechanicalskeleton of this robot is further edited based on the second templateof Fig. 6. We lengthened the torso and shortened the limbs to fitthis template into the input model (125 cm × 47 cm × 46 cm). Oursystem also produces plausible motion plans for this quad-robot. Failure case : While our system is stable in the generation ofslow walking motion, we find the generation of fast trot motion issensitive to the physical parameters. When the mass of the monsteris increased to two times, its influence on the mass center cannotbe balanced in the optimization and the COP constraint is violatedin the space-time optimization result, possibly due to its conflictbetween the foot contact constraint. Hence, the frame-by-frame UGUST 2019 11 Fig. 17. The monster-like robot takes two different input foot trajectories, and our system computes natural and physically correct motions for bothinputs.Fig. 18. Trotting motion for the monster-like robot. (a) Marked joint positions. The color of a dot indicates to which part of the horse skeleton itbelongs, and the positions are mapped to our monster skeleton joint angles using the space-time optimization with only kinematic constraints. (b)The designed trotting motion with our alternating algorithm. (c) The unstable motion simulated by the frame-by-frame optimization when the skindeformation is not considered in the space-time optimization.Fig. 19. The motion of a lizard-like robot. We edit the second template in Fig. 6 with SolidWorks scripts to create the design of its mechanicalskeleton. With user provided inputs, our system generates plausible motions of this robot. optimization will fail to produce a stable motion. Such situationmight be handled through the integration of foot plan samplingstep in [36]. ONCLUSION AND F UTURE W ORK In this paper, we have presented a fabrication-oriented motionplaning algorithm and detailed design/fabrication procedures forpersonalized skinned quad-robots. The physical constraints, suchas the equations of motion of the skinned robot and the motorconstraints, are integrated into the motion planning such thatthe resulting motion plan is physically and dynamically feasible. The condensation formulation allows us to conveniently establishthe nonlinear relationship between external forces and the targetkinematic parameters of the locomotion and to reach a QPCCformulation for the motion design. Our experiments show that thesystem is able to assist regular users to obtain natural and smoothmotions designed for skinned quad-robots.In the future, we want to explore a gait synthesis algorithmto generate the motion plan from high-level parameters, such asvelocity and turning angles. Combining captured gait data andoptimization with dynamic constraints has the potential to signifi-cantly reduce users’ labor efforts of creating such motion planning. UGUST 2019 12 Currently, the coupling between the FEM simulation of the softskin and the rigid body dynamics of the mechanical structure isnot fast enough for a closed-loop control of skinned quad-robots.We want to explore model reduction or homogenization techniquesto reduce the computational cost of FEM simulation and produceinteractive feedback to control skinned robots online. A PPENDIX AT HE E QUATIONS OF M OTION Lagrangian multibody dynamics : The multibody rigid skeletonof the robot is a kinematic tree of links connected by joints. Itsclassic Lagrangian mechanics formulation can be found in [84],where the DOFs of the skeleton are specified as the the generalizedcoordinate q of the joint angles and the root translation.The equation of motion for the articulated rigid skeleton canbe written as M r ( q ) ¨ q + D r ˙ q + f r ( q , ˙ q ) = g r , (20)where the subscript [ . ] r denotes variables for the skeleton rig.Detailed derivation of M r and f r in Eq. (20) can be found in theexcellent tutorial by Liu and Jain [85]. D r is the damping matrix.We refer to M r as the rigid mass matrix in order to differentiate itfrom the mass matrix of the soft skin. For a skeleton with K links, M r = K ∑ k = J (cid:62) k M ck J k , and M ck = (cid:20) m k · I 00 I ck (cid:21) , (21)where I is the identity matrix, I ck is the inertia tensor for the k th rigid link, and m k is the mass of the link. The Jacobian matrix J k = (cid:2) J (cid:62) vk , J (cid:62) ω k (cid:3) (cid:62) , where J vk = ∂ x k / ∂ q , relates the Cartesian coordinate x k of the link’s COM to the generalized coordinate q . Similarly, J ω k relates the angular velocity ωωω k to the generalized velocity ˙ q such that ωωω k = J ω k ˙ q . It is noteworthy that the rigid mass matrixis not constant because of the orientation-dependent inertia tensorand the Jacobian matrix.Non-inertia forces like Coriolis and centrifugal forces thatcouple the generalized coordinate are included in f r . The right-hand side of Eq. (20) is the generalized external force applied tothe skeleton, which includes the gravity force g and torques τττ fromthe actuating motors: g r = K ∑ k = J (cid:62) k (cid:20) g k τττ k (cid:21) = K ∑ k = (cid:104) J (cid:62) vk , J (cid:62) ω k (cid:105) (cid:20) g k τττ k (cid:21) . (22) FEM elastic simulation : The dynamics of the soft skin can alsobe formulated using Lagrangian mechanics, and we can obtain theequation of motion in a similar form: M d ¨ u + D d ˙ u + f d ( u ) = g d . (23)The subscript [ . ] d denotes variables for the deformable skin, and g d is the external force applied to the soft skin. We discretize thevolume of the skin by a tetrahedral mesh. The deformable massmatrix M d is constant and can be assembled using the standardFEM [86]. D d is the damping matrix, and f d denotes the internalelastic force, and it is the negative gradient of the strain energy Ψ such that f d = − ∇Ψ . The specific formulation of Ψ depends on thematerial model chosen. Typically, it is computed based on three isotropic invariants of the deformation gradient tensor F : I ( F ) = tr ( F (cid:62) F ) , I ( F ) = tr (cid:16) ( F (cid:62) F ) (cid:17) , I ( F ) = det ( F (cid:62) F ) . (24) For the robot with soft skin, its articulated skeleton motion willinduce large local compression, especially near the joint. Materialmodels like the StVK and co-rotational models that have beenwidely used in previous research [72], [74] become unstable underextreme compression. We therefore opted to use the Neo-Hookeanmaterial, whose strain energy density is defined as Ψ (cid:44) µ ( I − log ( I ) − ) + λ log ( I ) , (25)where µ and λ are Lam´e constants. The first Piola-Kirchhoff stresstensor P ∈ R × can be computed based on Eq. (25) using thechain rule: P = ∂ Ψ ∂ F = ∂ Ψ ∂ I · ∂ I ∂ F + ∂ Ψ ∂ I · ∂ I ∂ F = µ F − µ F −(cid:62) + λ log ( I ) F −(cid:62) . (26)The final formulation of f d is f d = − (cid:90) ∂ Ψ ∂ u d V = − (cid:90) ∂ Ψ ∂ F : ∂ F ∂ u d V = − (cid:90) (cid:18) P : ∂ F ∂ u (cid:19) (cid:62) d V . (27)Here, ∂ F / ∂ u ∈ R × × is a third-order tensor. For a tetrahedralelement with linear shape functions, ∂ F / ∂ u is constant and can beprecomputed and stored at each element. Implicit backward euler time integration : We discretize theequation of motions for both the rigid skeleton and deformablebody using the implicit backward Euler method to improve thestability of the simulation. Let ∆ q = q i + − q i and ∆ ˙ q = ˙ q i + − ˙ q i ,where the superscript [ · ] i is the frame index. Given the time interval ∆ t between frame i and i + 1, the velocity and acceleration of q at frame i + q i + = ∆ q / ∆ t and ¨ q i + = ∆ ˙ q / ∆ t = ∆ q / ∆ t − ˙ q i / ∆ t . The velocity and acceleration of u canbe derived similarly. Subsequently, Eqs. (20) and (23) can belinearized as (cid:0) M r ( q i ) + ∆ t D r + C ( q i , ˙ q i ) (cid:1) ∆ q = ∆ t (cid:0) g r − D r ˙ q i (cid:1) , (28)and (cid:18) M d + ∆ t D d + ∆ t ∂ f d ∂ u i (cid:19) ∆ u = ∆ t (cid:0) g d − f d ( u i ) − D d ˙ u i (cid:1) . (29)In Eq. (28), we compute M r ( q i ) and C ( q i , ˙ q i ) using state vari-ables at frame i . Therefore, Eq. (20) is only semi-implicitlydiscretized [87], [88]. Linear systems : The two-way coupled multibody-elastic systemin Eq.[2] in Sec. 4.1 of our paper is as follows: A r ∇ q C (cid:62) d ∇ u C (cid:62) ∇ q C ∇ u C ∆ q ∆ u λλλ = b r b d . (30)The matrices A r , A d in this equation are just the abbreviation ofsystem matrices in Eqs. (28) and (29), where A r = M r ( q i )+ ∆ t D r + ∆ t C ( q i , ˙ q i ) and A d = M d + ∆ t D d + ∆ t ∂ f d / ∂ u i . Analogically, wehave b r = ∆ t (cid:0) g r − D r ˙ q i − C ( q i , ˙ q i ) ˙ q i (cid:1) , and b d = ∆ t (cid:0) g d − f d ( u i ) − D d ˙ u i (cid:1) . A PPENDIX BS YSTEM M ATRIX C ONDENSATION The system matrix condensation in Sec.5.2 starts with eliminatingthe Lagrange multipliers λ in Eq. (30) (the same equation withEq.[2] in Sec. 4.1). We first expand the second line of this equation,which yields: A d ∆ u + ∇ u C (cid:62) λλλ = b d , or ∆ u = A − d (cid:0) b d − ∇ u C (cid:62) λλλ (cid:1) . (31) UGUST 2019 13 Similarly, expanding the first line in Eq. (30), we can produce anequation similar to Eq. (31) but for skeleton DOFs: ∆ q = A − r (cid:16) b r − ∇ q C (cid:62) λλλ (cid:17) . (32)Expanding the third line of Eq. (30) gives the linearized positionconstraint: ∇ q C ∆ q + ∇ u C ∆ u = . (33)Substituting both Eqs. (31) and (32) into Eq. (33) yields: λλλ = A − C (cid:0) ∇ q C A − r b r + ∇ u C A − d b d (cid:1) , (34)where A C = ∇ q C A − r ∇ q C (cid:62) + ∇ u C A − d ∇ u C (cid:62) . By substitutingEq. (34) back into Eqs. (31) and (32), we obtain the condensedformulas in Eq.[10] in our paper. R EFERENCES [1] C. Leger, “Automated synthesis and optimisation of robot configurations:An evolutionary approach,” Ph.D. dissertation, The Robotics Institute,Carnegie Mellon University, Pittsbugh, PA 15213, USA, 1999, cMU-RI-TR-99-43.[2] R. Desai, Y. Yuan, and S. Coros, “Computational abstractions forinteractive design of robotic devices,” in , May 2017, pp. 1196–1203.[3] V. Megaro, B. Thomaszewski, M. Nitti, O. Hilliges, M. Gross, andS. Coros, “Interactive design of 3d-printable robotic creatures,” ACMTrans. Graph. , vol. 34, no. 6, pp. 216:1–216:9, Oct. 2015.[4] P. Gallagher and M. Maclachlan, “Adjustment to an artificial limb: aqualitative perspective,” Journal of health psychology , vol. 6, no. 1, pp.85–100, 2001.[5] C. Majidi, “Soft robotics: a perspective - current trends and prospects forthe future,” Soft Robotics , vol. 1, no. 1, pp. 5–11, 2014.[6] Z. Wang, G. Hang, J. Li, Y. Wang, and K. Xiao, “A micro-robot fishwith embedded sma wire actuated flexible biomimetic fin,” Sensors andActuators A: Physical , vol. 144, no. 2, pp. 354–360, 2008.[7] J. Tan, G. Turk, and C. K. Liu, “Soft body locomotion,” ACM Trans.Graph. , vol. 31, no. 4, pp. 26:1–26:11, Jul. 2012.[8] X. Chen, H. Li, C.-W. Fu, H. Zhang, D. Cohen-Or, and B. Chen, “3dfabrication with universal building blocks and pyramidal shells,” in SIGGRAPH Asia 2018 Technical Papers , ser. SIGGRAPH Asia ’18,2018, pp. 189:1–189:15.[9] C. Dai, C. C. L. Wang, C. Wu, S. Lefebvre, G. Fang, and Y.-J. Liu,“Support-free volume printing by multi-axis motion,” ACM Trans. Graph. ,vol. 37, no. 4, pp. 134:1–134:14, Jul. 2018.[10] Y. Lan, Y. Dong, F. Pellacini, and X. Tong, “Bi-scale appearancefabrication,” ACM Trans. Graph. , vol. 32, no. 4, pp. 145:1–145:12, Jul.2013.[11] D. Chen, D. I. W. Levin, P. Didyk, P. Sitthi-Amorn, and W. Matusik,“Spec2fab: A reducer-tuner model for translating specifications to 3dprints,” ACM Trans. Graph. , vol. 32, no. 4, pp. 135:1–135:10, Jul. 2013.[12] A. Brunton, C. A. Arikan, T. M. Tanksale, and P. Urban, “3d printingspatially varying color and translucency,” ACM Trans. Graph. , vol. 37,no. 4, pp. 157:1–157:13, Jul. 2018.[13] K. Sakurai, Y. Dobashi, K. Iwasaki, and T. Nishita, “Fabricating reflectorsfor displaying multiple images,” ACM Trans. Graph. , vol. 37, no. 4, pp.158:1–158:10, Jul. 2018.[14] B. Bickel, M. B¨acher, M. A. Otaduy, H. R. Lee, H. Pfister, M. Gross, andW. Matusik, “Design and fabrication of materials with desired deformationbehavior,” ACM Trans. Graph. , vol. 29, no. 4, pp. 63:1–63:10, Jul. 2010.[15] M. Skouras, B. Thomaszewski, S. Coros, B. Bickel, and M. Gross,“Computational design of actuated deformable characters,” ACM Trans.Graph. , vol. 32, no. 4, pp. 82:1–82:10, Jul. 2013.[16] J. Panetta, Q. Zhou, L. Malomo, N. Pietroni, P. Cignoni, and D. Zorin,“Elastic textures for additive fabrication,” ACM Trans. Graph. , vol. 34,no. 4, pp. 135:1–135:12, Jul. 2015.[17] D. Chen, D. I. W. Levin, W. Matusik, and D. M. Kaufman, “Dynamics-aware numerical coarsening for fabrication design,” ACM Transactionson Graphics , vol. 36, no. 4, pp. 1–15, jul 2017.[18] M. Lau, A. Ohgawara, J. Mitani, and T. Igarashi, “Converting 3d furnituremodels to fabricatable parts and connectors,” ACM Trans. Graph. , vol. 30,no. 4, pp. 85:1–85:6, Jul. 2011. [19] J. Cal`ı, D. A. Calian, C. Amati, R. Kleinberger, A. Steed, J. Kautz, andT. Weyrich, “3d-printing of non-assembly, articulated models,” ACMTrans. Graph. , vol. 31, no. 6, pp. 130:1–130:8, Nov. 2012.[20] M. B¨acher, B. Bickel, D. L. James, and H. Pfister, “Fabricating articulatedcharacters from skinned meshes,” ACM Trans. Graph. , vol. 31, no. 4, pp.47:1–47:9, Jul. 2012.[21] S. Coros, B. Thomaszewski, G. Noris, S. Sueda, M. Forberg, R. W.Sumner, W. Matusik, and B. Bickel, “Computational design of mechanicalcharacters,” ACM Trans. Graph. , vol. 32, no. 4, pp. 83:1–83:12, Jul. 2013.[22] D. Ceylan, W. Li, N. J. Mitra, M. Agrawala, and M. Pauly, “Designingand fabricating mechanical automata from mocap sequences,” ACM Trans.Graph. , vol. 32, no. 6, pp. 186:1–186:11, Nov. 2013.[23] R. Zhang, T. Auzinger, D. Ceylan, W. Li, and B. Bickel, “Functionality-aware retargeting of mechanisms to 3d shapes,” ACM Trans. Graph. ,vol. 36, no. 4, pp. 81:1–81:13, Jul. 2017.[24] V. Megaro, J. Zehnder, M. Bcher, S. Coros, M. Gross, andB. Thomaszewski, “A computational design tool for compliant mech-anisms,” ACM Transactions on Graphics , vol. 36, no. 4, pp. 1–12, jul2017.[25] M. Geilinger, R. Poranne, R. Desai, B. Thomaszewski, and S. Coros,“Skaterbots: Optimization-based design and motion synthesis for roboticcreatures with legs and wheels,” ACM Trans. Graph. , vol. 37, no. 4, pp.160:1–160:12, Jul. 2018.[26] B. Bickel, P. Kaufmann, M. Skouras, B. Thomaszewski, D. Bradley,T. Beeler, P. Jackson, S. Marschner, W. Matusik, and M. Gross, “Physicalface cloning,” ACM Trans. Graph. , vol. 31, no. 4, pp. 118:1–118:10, Jul.2012.[27] J. M. Bern, K.-H. Chang, and S. Coros, “Interactive design of animatedplushies,” ACM Transactions on Graphics , vol. 36, no. 4, pp. 1–11, jul2017.[28] L.-K. Ma, Y. Zhang, Y. Liu, K. Zhou, and X. Tong, “Computational designand fabrication of soft pneumatic objects with desired deformations,” ACMTransactions on Graphics , vol. 36, no. 6, pp. 1–12, nov 2017.[29] A. Witkin and M. Kass, “Spacetime constraints,” SIGGRAPH Comput.Graph. , vol. 22, no. 4, pp. 159–168, Jun. 1988.[30] J. V. Albro, G. A. Sohl, J. E. Bobrow, and F. C. Park, “On the computationof optimal high-dives,” in IEEE International Conference on Roboticsand Automation. , vol. 4, 2000, pp. 3958–3963.[31] M. F. Cohen, “Interactive spacetime control for animation,” SIGGRAPHComput. Graph. ACM Trans. Graph. , vol. 22, no. 3, pp. 417–426, Jul.2003.[34] A. Safonova, J. K. Hodgins, and N. S. Pollard, “Synthesizing physicallyrealistic human motion in low-dimensional, behavior-specific spaces,” ACM Trans. Graph. , vol. 23, no. 3, pp. 514–521, Aug. 2004.[35] A. Sulejmanpaˇsi´c and J. Popovi´c, “Adaptation of performed ballisticmotion,” ACM Trans. Graph. , vol. 24, no. 1, pp. 165–179, Jan. 2005.[36] K. Wampler and Z. Popovi´c, “Optimal gait and form for animallocomotion,” ACM Trans. Graph. , vol. 28, no. 3, pp. 60:1–60:8, Jul.2009.[37] X. Wei, J. Min, and J. Chai, “Physically valid statistical models for humanmotion generation,” ACM Transactions on Graphics , vol. 30, no. 3, pp.1–10, may 2011.[38] K. Wampler, Z. Popovi´c, and J. Popovi´c, “Generalizing locomotion styleto new animals with inverse optimal regression,” ACM Trans. Graph. ,vol. 33, no. 4, pp. 49:1–49:11, Jul. 2014.[39] Z. Popovi´c and A. Witkin, “Physically based motion transformation,” in Proceedings of the 26th Annual Conference on Computer Graphics andInteractive Techniques , ser. SIGGRAPH ’99, 1999, pp. 11–20.[40] M. H. Raibert and J. K. Hodgins, “Animation of dynamic leggedlocomotion,” SIGGRAPH Comput. Graph. , vol. 25, no. 4, pp. 349–358,Jul. 1991.[41] J. K. Hodgins, W. L. Wooten, D. C. Brogan, and J. F. O’Brien, “Animatinghuman athletics,” in Proceedings of the 22Nd Annual Conference onComputer Graphics and Interactive Techniques , ser. SIGGRAPH ’95.ACM, 1995, pp. 71–78.[42] S. Coros, A. Karpathy, B. Jones, L. Reveret, and M. van de Panne,“Locomotion skills for simulated quadrupeds,” ACM Trans. Graph. , vol. 30,no. 4, pp. 59:1–59:12, Jul. 2011.[43] M. Hirofumi and I. Shimoyama, “Dynamic walk of a biped,” TheInternational Journal of Robotics Research , vol. 3, no. 2, pp. 60–74,1994. UGUST 2019 14 [44] J. K. Hodgins and N. S. Pollard, “Adapting simulated behaviors for newcharacters,” in Proceedings of the 24th Annual Conference on ComputerGraphics and Interactive Techniques , ser. SIGGRAPH ’97, 1997, pp.153–162.[45] A. D. Kuo, “Stabilization of lateral motion in passive dynamic walking,” The International Journal of Robotics Research , vol. 18, no. 9, pp. 917–930, 1999.[46] S. Kudoh, T. Komura, and K. Ikeuchi, “Stepping motion for a human-likecharacter to maintain balance against large perturbations,” in Proceedings2006 IEEE International Conference on Robotics and Automation, 2006.ICRA 2006. , 2006, pp. 2661–2666.[47] T. Niwa, S. Inagaki, and T. Suzuki, “Locomotion control of multi-leggedrobot based on follow-the-contact-point gait,” in , 2009,pp. 2247–2253.[48] Y. Li, B. Li, J. Ruan, and X. Rong, “Research of mammal bionic quadrupedrobots: A review,” in , Sept 2011, pp. 166–171.[49] K. Yin, K. Loken, and M. van de Panne, “Simbicon: Simple bipedlocomotion control,” ACM Trans. Graph. , vol. 26, no. 3, Jul. 2007.[50] L. Liu, K. Yin, M. van de Panne, T. Shao, and W. Xu, “Sampling-basedcontact-rich motion control,” ACM Trans. Graph. , vol. 29, no. 4, pp.128:1–128:10, Jul. 2010.[51] S. Ha, Y. Ye, and C. K. Liu, “Falling and landing motion control forcharacter animation,” ACM Trans. Graph. , vol. 31, no. 6, pp. 155:1–155:9,Nov. 2012.[52] L. Liu, K. Yin, M. van de Panne, and B. Guo, “Terrain runner: Control,parameterization, composition, and planning for highly dynamic motions,” ACM Trans. Graph. , vol. 31, no. 6, pp. 154:1–154:10, Nov. 2012.[53] M. da Silva, Y. Abe, and J. Popovi´c, “Interactive simulation of stylizedhuman locomotion,” ACM Trans. Graph. , vol. 27, no. 3, pp. 82:1–82:10,Aug. 2008.[54] U. Muico, Y. Lee, J. Popovi´c, and Z. Popovi´c, “Contact-aware nonlinearcontrol of dynamic characters,” ACM Trans. Graph. , vol. 28, no. 3, pp.81:1–81:9, Jul. 2009.[55] S. Jain and C. K. Liu, “Controlling physics-based characters using softcontacts,” ACM Trans. Graph. , vol. 30, no. 6, pp. 163:1–163:10, Dec.2011.[56] J. Kim and N. S. Pollard, “Fast simulation of skeleton-driven deformablebody characters,” ACM Trans. Graph. , vol. 30, no. 5, pp. 121:1–121:19,Oct. 2011.[57] T. Shinar, C. Schroeder, and R. Fedkiw, “Two-way coupling ofrigid and deformable bodies,” in Proceedings of the 2008 ACM SIG-GRAPH/Eurographics Symposium on Computer Animation , ser. SCA ’08,2008, pp. 95–103.[58] C. Duriez, “Control of elastic soft robots based on real-time finiteelement method,” in , May 2013, pp. 3982–3987.[59] F. Largilliere, V. Verona, E. Coevoet, M. Sanz-Lopez, J. Dequidt, andC. Duriez, “Real-time control of soft-robots using asynchronous finiteelement modeling,” in , May 2015, pp. 2550–2555.[60] C. Duriez and T. Bieze, “Soft robot modeling, simulation and control inreal-time,” in Soft Robotics: Trends, Applications and Challenges , 2017,pp. 103–109.[61] E. Coevoet, T. Morales-Bieze, F. Largilliere, Z. Zhang, M. Thieffry,M. Sanz-Lopez, B. Carrez, D. Marchal, O. Goury, J. Dequidt, andC. Duriez, “Software toolkit for modeling, simulation, and control ofsoft robots,” Advanced Robotics , vol. 31, no. 22, pp. 1208–1224, 2017.[62] A. K. Z. Z. R. M. Thor Morales Bieze, Frederick Largilliere and C. Duriez,“Finite element method-based kinematics and closed-loop control of soft,continuum manipulators,” Soft Robotics , vol. 5, no. 3, 2018.[63] D. Terzopoulos, J. Platt, A. Barr, and K. Fleischer, “Elastically deformablemodels,” SIGGRAPH Comput. Graph. , vol. 21, no. 4, pp. 205–214, 1987.[64] G. Irving, J. Teran, and R. Fedkiw, “Tetrahedral and hexahedral invertiblefinite elements,” Graphical Models , vol. 68, no. 2, pp. 66–89, 2006.[65] G. Hirota, S. Fisher, and M. Lin, “Simulation of non-penetrating elasticbodies using distance fields,” Chapel Hill, NC, USA, Tech. Rep., 2000.[66] M. M¨uller, J. Dorsey, L. McMillan, R. Jagnow, and B. Cutler, “Sta-ble real-time deformations,” in Proceedings of the 2002 ACM SIG-GRAPH/Eurographics symposium on Computer animation . ACM, 2002,pp. 49–54.[67] A. Nealen, M. M¨uller, R. Keiser, E. Boxerman, and M. Carlson,“Physically based deformable models in computer graphics,” vol. 25,no. 4, pp. 809–836, 2006.[68] J. Barbiˇc, M. da Silva, and J. Popovi´c, “Deformable object animationusing reduced optimal control,” ACM Trans. Graph. , vol. 28, no. 3, pp.53:1–53:9, Jul. 2009. [69] J. Barbiˇc, F. Sin, and E. Grinspun, “Interactive editing of deformablesimulations,” ACM Trans. Graph. , vol. 31, no. 4, pp. 70:1–70:8, Jul. 2012.[70] S. Li, J. Huang, F. de Goes, X. Jin, H. Bao, and M. Desbrun, “Space-timeediting of elastic motion through material optimization and reduction,” ACM Trans. Graph. , vol. 33, no. 4, pp. 108:1–108:10, Jul. 2014.[71] Z. Pan and D. Manocha, “Active animations of reduced deformablemodels with environment interactions,” ACM Trans. Graph. ,vol. 37, no. 3, pp. 36:1–36:17, Aug. 2018. [Online]. Available:http://doi.acm.org/10.1145/3197565[72] S. Capell, S. Green, B. Curless, T. Duchamp, and Z. Popovi´c, “Interactiveskeleton-driven dynamic deformations,” in ACM Transactions on Graphics(TOG) , vol. 21, no. 3. ACM, 2002, pp. 586–593.[73] S.-H. Lee, E. Sifakis, and D. Terzopoulos, “Comprehensive biomechanicalmodeling and simulation of the upper body,” ACM Transactions onGraphics (TOG) , vol. 28, no. 4, p. 99, 2009.[74] L. Liu, K. Yin, B. Wang, and B. Guo, “Simulation and control of skeleton-driven soft body characters,” ACM Transactions on Graphics (TOG) ,vol. 32, no. 6, p. 215, 2013.[75] R. Ogden, “Large deformation isotropic elasticity-on the correlationof theory and experiment for incompressible rubberlike solids,” in Proceedings of the Royal Society of London A: Mathematical, Physicaland Engineering Sciences , vol. 326, no. 1567. The Royal Society, 1972,pp. 565–584.[76] M. Pozzi, E. Miguel, R. Deimel, M. Malvezzi, B. Bickel, O. Brock, andD. Prattichizzo, “Efficient fem-based simulation of soft robots modeled askinematic chains,” in , May 2018, pp. 1–8.[77] M. Moore and J. Wilhelms, “Collision detection and response for computeranimation,” in Proceedings of the 15th Annual Conference on ComputerGraphics and Interactive Techniques , ser. SIGGRAPH ’88, 1988, pp.289–298.[78] R. Bridson, R. Fedkiw, and J. Anderson, “Robust treatment of collisions,contact and friction for cloth animation,” ACM Trans. Graph. , vol. 21,no. 3, pp. 594–603, Jul. 2002.[79] M. Bro-Nielsen and S. Cotin, “Real-time volumetric deformable modelsfor surgery simulation using finite elements and condensation,” ComputerGraphics Forum , vol. 15, no. 3, pp. 57–66, 1996.[80] Y. Teng, M. Meyer, T. DeRose, and T. Kim, “Subspace condensation:Full space adaptivity for subspace deformations,” ACM Trans. Graph. ,vol. 34, no. 4, pp. 76:1–76:9, Jul. 2015. [Online]. Available:http://doi.acm.org/10.1145/2766904[81] M. Anitescu and F. A. Potra, “Formulating dynamic multi-rigid-body con-tact problems with friction as solvable linear complementarity problems,” Nonlinear Dynamics , vol. 14, no. 3, pp. 231–247, 1997.[82] O. K.-C. Au, C.-L. Tai, H.-K. Chu, D. Cohen-Or, and T.-Y. Lee, “Skeletonextraction by mesh contraction,” ACM Trans. Graph. , vol. 27, no. 3, pp.44:1–44:10, Aug. 2008.[83] T.-C. Huang, Y.-J. Huang, and W.-C. Lin, “Real-time horse gait synthesis,” Computer Animation and Virtual Worlds , vol. 24, no. 2, pp. 87–95, 2013.[84] A. A. Shabana, Dynamics of multibody systems . Cambridge universitypress, 2013.[85] C. K. Liu and S. Jain, “A quick tutorial on multibody dynamics,” Onlinetutorial, June , p. 7, 2012.[86] K.-J. Bathe, Finite element method . Wiley Online Library, 2008.[87] D. Baraff and A. Witkin, “Large steps in cloth simulation,” in Proceedingsof the 25th Annual Conference on Computer Graphics and InteractiveTechniques , ser. SIGGRAPH ’98. New York, NY, USA: ACM, 1998, pp.43–54.[88] P. Song, J.-S. Pang, and V. Kumar, “A semi-implicit time-stepping modelfor frictional compliant contact problems,” International Journal forNumerical Methods in Engineering , vol. 60, no. 13, pp. 2231–2261, 2004. Xudong Feng Xudong Feng is a Ph.D. candidatein the State Key Lab of CAD & CG, Colledgeof Computer Science, Zhejiang University. Hereceived his bachelor degree in Tianjin University.His research interests include physical basedsimulation and animation, digital fabrication andreinforcement learning. UGUST 2019 15 Jiafeng Liu is currently working toward the Mas-ter degree in computer engineering with theZhejiang University. He received the bachelor’sdegree from the Hefei University of Technology in2018. He His research interests include characteranimation, robotics locomotion, machine learningand related topics. Huamin Wang received the BEng degree fromZhejiang University, the MS degree from StanfordUniversity, and the PhD degree in computerscience from the Gerogia Institute of Technology,in 2002, 2004 and 2009. Hw is an associateprofessor in the Department of Computer Sci-ence and Engineering, the Ohio State University.Before joining Ohio State University, he was apostdoctoral researcher in the Department ofElectrical Engineering and Computer Sciences,the University of California, Berkeley. He is amember of the IEEE. Yin Yang received the PhD degree in computerscience from the University of Texas at Dallas,in 2013. He is an assistant professor in theDepartment of Electrical Computer Engineering,the University of New Mexico, Albuquerque. Hisresearch interests include physics-based anima-tion/ simulation and related applications, scientificvisualization, and medical imaging analysis. Heis a member of the IEEE. Hujun Bao is a Chenkong professor in State KeyLab of CAD&CG, College of Computer Science atZhejiang University. He received PhD degree inapplied mathematics from Zhejiang university in1993. His main research interests are computergraphics and computer vision, including real-timerendering technique, geometry computing, virtualreality, and 3D reconstruction, and has publishedmore than 100 papers on prestigious academicjournals and international conferences. Bernd Bickel is an Assistant Professor, headingthe Computer Graphics and Digital Fabricationgroup at IST Austria. He received Master andPh.D degree in Computer Science from ETHZurich. His research interests includeds computergraphics and its overlap into robotics, computervision, biomechanics, material science, and digi-tal fabrication. He receives SIGGRPAH significantresearcher award in 2017.