Physical Cyclic Animations

We address the problem of synthesizing physical animations that can loop seamlessly. We formulate a variational approach by deriving a physical law in a periodic time domain. The trajectory of the animation is represented as a parametric closed curve, and the physical law corresponds to minimizing the bending energy of the curve. Compared to traditional keyframe animation approaches, our formulation is constraint-free, which allows us to apply a standard Gauss--Newton solver. We further propose a fast projection method to efficiently generate an initial guess close to the desired animation. Our method can handle a variety of physical cyclic animations, including clothes, soft bodies with collisions, and N-body systems.


INTRODUCTION
Cyclic animations are short animations that can be played in loops.They have existed since the early history of animations in the form of phenakistoscopes (Fig. 1).Looping animations continue to draw great interest in modern days for video games and movies. 1 They also play a central role in keyframe-controlled animations [Barbič et al. 2009], and designing characters' and animatronics' cyclic motion such as walking [Nunes et al. 2012;Mordatch et al. 2013;Starke et al. 2022].
We aim to produce cyclic animations of physical systems.The animation should locally approximate the physical laws of motion, and globally loop back seamlessly (Fig. 2).In particular, we aim to reduce the number of user inputs, such as keyframes for guiding the animation [Barbič et al. 2009].Our problem setting poses new challenges: while fewer keyframe fidelity constraints leave more room for chaotic physical systems to develop a larger variety of interesting motions, it produces worse control gradients evaluated by the common adjoint methods.Fig. 2. Flag post.Given a physical system with an initial condition (top), our method synthesizes a physical animation (bottom) that can loop back seamlessly, such that the end frame would be the same as the initial frame.Here we show a cyclic animation of a flag under the wind.
We describe a variational approach to generate physical cyclic animations.The cyclic simulation is represented as a parameterized closed curve (loop) in the configuration space of the physical system (Fig. 3).The physical law assigns to each configuration point an associated acceleration of the trajectory.A trajectory passing through a point with a different acceleration is similar to a bent curve which gives rise to a cost of bending energy.Using this interpretation of a residual force minimization problem, we model our cyclic physical simulation as a closed curve of minimal total bending energy.Users can specify initial conditions (positions/velocities at the starting frame) either as soft or hard constraints of specific points of the loop.
Our closed curve parameterization guarantees a cyclic animation.Moreover, the resulting minimization problem is essentially unconstrained, in which case many optimization methods directly apply.We present a fast projection method for an initial guess followed by a Gauss-Newton method to obtain the optimal solution.Our optimization can often overcome the challenges caused by the non-linear constraints when applying the classical shooting method to a chaotic system (Fig. 3).
We show results on a variety of physical systems with a moderately large degrees of freedom , including thin shell models, hyperelastic finite element models with contact, and N-body simulation (see Section 5).
Comparing with the conventional key-frame animation control methods, we highlight: • We build the cyclic constraint in the topology of trajectory, leaving an essentially unconstrained optimization problem.Solve for an initial guess Gauss-Newton solve Fig. 3.For a chaotic dynamical system, a direct gradient descent method for the optimal control problem (9) often diverges due to the amplification in the sensitivity propagation (b).Our method (c) overcomes the problem by treating the time domain as a closed loop and avoiding using gradient descent for the initial guess.
• In the traditional optimal control approach, the gradient descent can be unstable or inefficient since the trajectory can be sensitive to control force in a general dynamical system.We replace the gradient descent procedure with a more efficient and robust fast projection method.
• The heuristic fast projection is backed by a Gauss-Newton solve that seeks a local minimum.
The method works even for chaotic dynamical systems where traditional gradient-based optimal control fails.

RELATED WORK
Some current practice for generating physical looping animations is to interpolate between states of the end frame and the initial frame, or to resort to non-physical periodic motions such as a sine wave (Footnote 1).These approaches guarantee cyclic animations, at the cost of producing less physical motion.Some approaches solve partial differential equations with periodic boundary conditions using spectral methods [Castro et al. 2010].
Spacetime constraints.Our method belongs to the methods of spacetime constraints [Witkin and Kass 1988].These methods specify the initial and end conditions (the keyframes) and find the most physically accurate trajectories that satisfy the conditions by applying control forces.In its general form, spacetime constraints involve non-convex optimization with non-linear constraints.Prior work has demonstrated the use of spacetime constraints to generate periodic motions [Barbič et al. 2009;Wampler and Popović 2009;Nunes et al. 2012;Schulz et al. 2014].
While spacetime constraints can be solved using techniques such as gradient descent and sequential quadratic programming, meeting the non-linear constraints is hard.The optimization becomes especially challenging when the keyframes are distant apart, leaving long time intervals that can support chaotic dynamics.
To address the challenges in optimization and provide user control, previous work employed several strategies.First, multiple keyframes are often used to further constrain the animation.Next, the non-linear optimization is often projected to a lower-dimensional subspace [Barbič et al. 2009] or converted into a linear system [Huang et al. 2011;Barbič et al. 2012;Hildebrandt et al. 2012;Li et al. 2014;Schulz et al. 2014].Alternatively, some work applies derivative-free optimization [Wampler and Popović 2009;Nunes et al. 2012].
Our method targets degrees of freedom higher than usual optimal control approaches (e.g., we optimize for control forces on each vertex of a cloth), while focusing on the case where the only constraints are that the animation should loop back and be as physical as possible.Our method does not require additional keyframes (other than the initial condition), subspace projection, or linearization, that could sacrifice the physics.This is enabled by our constraint-free reformulation of the spacetime constraints problem and our Gauss-Newton optimization, whereas previous work still needs to meet the constraints even with the aforementioned remedies.On the other hand, we leave having more user controls as future work.
Shooting methods and optimal control.A common way of solving spacetime constraints in graphics is the shooting method [Betts 1998]: given a forward simulation, the constrained optimization is solved using the derivatives with respect to the control forces [Popović et al. 2003;Treuille et al. 2003;Wojtan et al. 2006].The derivatives are often computed by the adjoint method [Johnson 2012], which shares similar computational efficiency to backpropagation or reverse-mode automatic differentiation.Some recent differentiable physics work also belong to this category [de Avila Belbute-Peres et al. 2018;Hu et al. 2020;Du et al. 2021], as well as work on trajectory optimization [Zimmermann et al. 2019], physical design [Zehnder et al. 2021], and fluid control [McNamara et al. 2004;Pan et al. 2013;Tang et al. 2021].
Cyclic animation can also be solved using the shooting method, by placing a constraint on the ending state.Prior work that demonstrates periodic motions using spacetime constraints often solves the optimization using the shooting method [Barbič et al. 2009;Schulz et al. 2014].However, meeting the non-linear constraints through optimizing control forces is challenging.Another class of optimal control methods, direct collocation [Hargraves and Paris 1987;Posa et al. 2014], optimizes for both the control forces and trajectory, but they still need to handle non-linear constraints.We show that shooting methods, when used for generating cyclic animations, can violate the constraints and produce non-smooth "jumps" in the animation.Our method optimizes only the trajectory, making it trivial to satisfy the constraints, and can be used in conjunction with the shooting methods to improve their results.
Image/video-based methods.Some methods directly process videos to generate infinite loops [Schödl et al. 2000;Schödl and Essa 2002;Liao et al. 2013;He et al. 2017;Härkönen et al. 2022].Among this class of methods, some use physically-inspired periodic motion to produce cyclic videos [Chuang et al. 2005;Davis et al. 2015].We instead directly operate on the 3D physical states.
Data-driven motion synthesis.Some recent work applies data-driven techniques to synthesize a variety of (potentially periodic) motions from examples [Starke et al. 2022;Li et al. 2022].We focus on physics-based approaches, and leave the combination between our work and data-driven methods as future work.

FORMULATION
Our goal is to formulate a physical law on a periodic time domain.We consider physical systems whose equations of motion take the form: where R  represents the physical configuration space,2 M = (  ) ∈ R × the mass, and F a known force law.Our problem of producing a cyclic animation takes a set of desired initial values q(0) = r 0 , q(0) = v 0 and a time period  > 0 as inputs.
To model temporal periodicity, we let the domain of time be the periodic interval [0, ), which we denote by the quotient space (2) We consider only those continuous paths q : S 1  → R  in the configuration space parameterized by this cyclic time domain.Define the feasible set as the vector space of parameterized closed curves (Fig. 3) For each closed curve q ∈ V we define the energy where the integrand  is a "bending energy" for the curve that measures its deviation from satisfying the physical system (1) where |u| 2 Such a least-squares formulation follows Gauss' principle of least constraint [Gauß 1829;Hildebrandt et al. 2012].

Optimization problems
Our main optimization problem is to find q ∈ V that solves Here, ,  ′ are parameters that control how much we desire the fidelity for the initial values.
When ,  ′ → 0, these "soft" initial conditions become hard linear constraints.In this limit, the feasible set reduces to the affine subspace over which we solve minimize Fig. 4. Starting with a forward simulation of a 3-body system (a), we compare solving the optimal control problem (9) with non-linear constraints, using soft constraints (b), and projected gradient descent (c).The figure visualizes the trajectories.Projected gradient descent only converged for 9% of the initial conditions in our randomized trials: here we show the rare case where it converges.Our constraint-free formulation allows us to generate better seamless and physical animations (d) compared to the traditional methods based on solving non-linear constraints.

Relation to optimal control
With a substitution u = M q − F(q), our optimization is equivalent to an optimal control problem: The optimizer u is the minimal control force added to the system (1) to guide q to meet the constraints at  =  .Remark 1.Although the optimal control problem (9) is equivalent to the unconstrained problem (8) by variable substitution, the non-linear physical constraint M q = F(q) + u makes satisfying the cyclic constraints much more difficult.Constraint violation is not tolerable for seamless cyclic animation.Constrained optimization, such as Lagrange multipliers or soft constraints with gradient descent, works well for small time intervals with monotonic motion.Such control becomes much more difficult for animation over a long time interval that supports complex physical motions.In practice, for the optimal control approach to work robustly, one adds multiple keyframes throughout the animation interval (rather than just an initial condition) [Barbič et al. 2009]; however, each additional keyframe constraint can sacrifice the physics.Our method builds cyclicity into the problem and removes the cyclic constraint from optimization.
Our experiments (Fig. 4c) show that solving (9) with projected gradient descent or Lagrange multipliers usually leads to divergent iterations with increased loss and violated constraints [Platt and Barr 1988].It is possible to relax the constraints q(0) = q( ) and q(0) = q( ) into soft constraints [Zehnder et al. 2021;Du et al. 2021], but this can lead to a sudden jump between the end and initial states (Fig. 4b).Optimizing soft constraints also requires extensive parameter tuning.

METHOD
We develop an algorithm to solve the variational problem (6).We discretize the periodic time domain and formulate the discrete counterpart of the optimization (6) (Section 4.1), which we solve iteratively using the Gauss-Newton method (Section 4.2) with an efficient initial guess (Section 4.3).
We use index notation to specify the components of the variable.An element  ()   is the scalar value of the -th component of the physical configuration at time step  at optimization iteration .We use boldface q  = ( 0 ,  1 , • • • ,  ,−1 ) ∈ R  to denote the physical state at time step .

Discretization
We discretize the temporal domain (2) uniformly into  nodes We use arithmetic modulo  on the index  ∈ Z  , and denote time step size by ℎ =  / .Under this time-discretization the collection of closed curves (3) becomes a finite-dimensional vector space We define the operator that takes the second partial derivative with respect to time Using this operator we define the discrete version of the total loss function (4) into a functional  :  → R. We have two options of the discrete loss function depending on the time integration scheme: and ∈Z  ℎ q  , (q)  (14a) The difference between ( 13) and ( 14) corresponds to the two time integration schemes: backward Euler uses the force F(q +1 ) for the acceleration (q)  , where as symplectic Euler uses F(q  ). 3 With the discrete energy ((13) or ( 14)), our main optimization problem (6) becomes We add ℎ 3 factor in the denominator to balance the units.The desired initial positions r 0 , r 1 and the penalty parameters  0 ,  1 are prescribed by the user.In our discretization, the velocity q(0) is discretized using finite difference q 1 −q 0 ℎ .In the following, we stick with the symplectic energy ( 14).Converting to backward Euler is straightforward.

Optimization with Gauss-Newton method
We seek the minimum of (15) iteratively by the Gauss-Newton method, which is suitable for our loss function as a sum of squares 3 The choices for the arguments of ensure that the absolute zero energy state (q  +1 , (q)  ) = 0 (resp.(q  , (q is equivalent to the commonly employed backward (resp.symplectic) Euler time integration [Baraff and Witkin 1998] backward Proc.ACM Comput.Graph.Interact.Tech., Vol. 6, No. 2, Article .Publication date: August 2023.
Equation ( 16) can be written as a quadratic form: where B = ℎ and g :  → R 2 by 4.2.1 Gauss-Newton method.A Gauss-Newton method is a quasi-Newton method that uses a specific approximation of the Hessian of the loss function.The gradient of our loss function is given by where J = u q and K = g q are the Jacobians of ( 18) and ( 19) respectively (see Section 4.2.2 for the explicit formulas).In a Gauss-Newton method, one takes the semi-positive definite approximation for the Hessian This yields the following quasi-Newton algorithm.
Algorithm 1: Main algorithm input: Initial guess q (0) from Algorithm 2 1 while |q (+1) − q () | <  do 2 Evaluate u, g, J = u q and K = g q at q () ⊲ (18), ( 19) q (+1) ← q () + q ⊲  ∈ (0, 1) given by line search.7 end output : q () in the last iteration One can use standard methods to solve the linear system in Step 5. Our current implementation uses a Cholesky decomposition.4.2.2Jacobians of the residual forces.The Jacobian J = u q of the residual force ( 18) is given by for all variation q ∈  , where is the negative gradient of the force law evaluated at configuration q  ∈ R  .If F is a conservative force F = −∇ arising from a potential energy  : R  → R, then H  = Hess  | q  .The Jacobian K = g q of ( 19) is given by Kq = q0 q1 . Proc

Hard position constraints.
When  0 ,  1 → 0 or when there are boundary conditions in the physical system, the values of   are prescribed and fixed during optimization for a subset of ( , ) ∈ Z  ×Z  .In the presence of such hard position constraints, we simply remove the prescribed   's from the set of unknown variables.This procedure reduces the size of the matrix Hess , and, in the case of  0 ,  1 = 0, removes the stiff term K ⊺ CK.

Initial guess with fast projection
For the initial guess for Algorithm 1, we construct a heuristic q (0) that is close to being cyclic and physical.To achieve this, we solve an approximate solution to the control problem (9).Define a functional q that maps a control force u to the solution to the initial value problem (forward simulation) With this substitution, we rewrite (9) as where We apply the fast projection method [Goldenthal et al. 2007;Dinev et al. 2018] to find a point on the constraint manifold { = 0} close to the origin as an approximate solution to (25).As illustrated in Fig. 5, the method starts with u (0) = 0 and iteratively updates with the orthogonal projection   ( ) onto the zero-hyperplane  ( ) of the linearized constraint functional at u ( ) : The stepping u = u (+1) − u ( ) is perpendicular to the hyperplane  ( ) and therefore where  ∈ R is a scalar multiplier.The projection step ( 27) thus amounts to solving We can solve the linear system efficiently: observe (by Gauss elimination) that  =  /|∇ | 2 , and therefore we can evaluate u = −∇ = −  | ∇ | 2 ∇ explicitly.We evaluate ∇ by the standard adjoint method or reverse-mode automatic differentiation [Johnson 2012].The overall process of finding an initial guess only adds little computational overhead.We summarize our algorithm for finding an initial guess below.

⊲ forward simulation
Remark 2. The fast projection (30) is an approximation of the projected gradient descent method by assuming both the control force and the Lagrange multiplier are zero.This relaxation helps stabilizing the optimization when the system is less stiff.As a cautionary note, the iteration could diverge when the required control force is large and the approximation is no longer adequate.

RESULTS
We demonstrate that our algorithm generates physical cyclic animations for various physical systems, including N-body systems, cloth simulation under externally driven forces or periodic boundary conditions, and soft body simulation with contact.We then evaluate the effectiveness of fast projection as the initial guess.We also show one example where fast projection fails, but Gauss-Newton solver can still generate physical cyclic animation by taking forward simulation as an initial guess.
For cloth and soft body simulation, we adopt the finite element implementation from the C-IPC codebase [Li et al. 2021].We also implement penalty-based collision with external objects.We use CHOLMOD [Chen et al. 2008] for its linear solvers and OpenMP for parallel computing.All related experiments are run on a 2.3 GHz 8-Core Intel Core i7-11800H CPU with 32 GB memory.Apart from that, we implement the N-body example in Python and execute it on an 8-Core Apple M1 and 8 GB memory.We use symplectic integration in N-Body and implicit integration in cloth and soft body simulation for stability.We report the experiment details and optimization run time in Table 1.A time step of 10 ms is used in all of the experiments.Fig. 6.N-body.A 5-body system of equal mass interacting under mutual gravity.Physically, the system evolves in a chaotic manner (left).Our proposed method allows for synthesizing a realistic cyclic animation (right).The trajectories of forward simulation experience a jump (visualized in dashed lines), whereas our method loops back to the initial orbit.

N-body system
A number (N) of celestial bodies interacting under Newton's gravity often results in chaotic dynamics (for N ≥ 3) [ Barrow-Green 1997].Seeking time-periodic solutions in an N-body system has been a computationally challenging problem in astronomy [Li and Liao 2019].Despite the chaotic nature of the physical system, our algorithm can generate close-to-physical time-periodic solutions with relatively small computational cost.
As shown in Fig. 6, a physical simulation of an N-body problem by solving an initial value problem (left) is unlikely to be close to cyclic.Our method simulates the system while guaranteeing time periodicity.The trajectories of the 5 bodies form 5 interlinking rings.

Cloth
We demonstrate cyclic cloth simulations using a thin shell [Grinspun et al. 2003] with bending energy and Neo-Hookean membrane energy using the physical material parameters of 100% cotton.5.2.1 Flag Post.Fig. 2 shows an animation of a rectangular flag with two fixed corners.In addition to the elastic forces, we add a wind force proportional to the dot product between the surface normal and a constant wind vector. 4To get natural starting positions, we run forward simulation and use the configurations from the 60 th and 61 st frame to set r 0 and r 1 respectively.Forward simulations with this setup yield a physics-based animation of a waving flag (Fig. 2 left), but the animation generally does not repeat.Our initial guess (Section 4.3) finds a plausible solution, but still exhibits an unnatural jump at the end.Optimizing using our method from the initial guess finds a physically realistic cyclic animation that shows no visual artifact.
5.2.2 Ribbons.Fig. 7 shows an animation of narrow strips of thin shell approximating anisotropic elastic rods [Bergou et al. 2008].We attach three ribbons to a point handle which orbits in cyclic motion.Even with such a time-periodic external influence, the forward simulation fails to loop back.Our method is able to produce a natural motion of ribbons that also loops back seamlessly.

Soft body with collision
We apply our method to soft body simulation with collision.We adopt a penalty-based contact model with the floor and sphere [Macklin et al. 2020].

Bouncing Doughnut.
As shown in Fig. 8, we setup a doughnut-shaped soft body with a downward velocity at frame 0. In the forward simulation, the body bounces off to a different course after its collision with the ground and never loops back.With our initial guess, the body manages to loop back to its initial position.However, this approximated solution still suffers from a large virtual acceleration around the last frames, which yields artifacts in video playback.This is also visible as a discontinuity in the velocity curve (Fig. 9).Our optimal solution produces a natural Our method can be applied to soft body simulation with contact.Our initial guess generates an animation that almost loops back, at the cost of losing some momentum due to a large control force, making it not bounce as high.Our optimization produces a seamless loop, while preserving the momentum.
cyclic animation.The soft body bounces up to the correct height (frame 40), from which it falls back to the initial position, matching the velocity of the initial frame.Fig. 9 shows the vertical component of the center of mass plotted against time.The physical forward simulation loses the vertical momentum as it is transferred to an angular one upon collision.The fast projection is not able to restore this loss of vertical momentum.Our method restores the vertical momentum and successfully loops back in both position and velocity.5.3.2Wobbly bar.In Fig. 10, we show a scripted rigid ball moving through the scene in constant velocity and colliding with an elastic bar.The elastic bar is initialized from the rest shape.This is a challenging scenario with large deformation and a long time interval.The forward simulation fails to loop back, while our method generates natural and seamless animation.

Swinging cloth.
Here we compare fast projection with gradient descent in solving keyframe control problems and evaluate their role as initial guesses.We found that directly solving (9) using projected gradient descent usually diverges.Instead, previous spacetime constraints methods [Barbič et al. 2009;Wojtan et al. 2006]  A piece of cloth swings back and forth from a horizontal initial state.The dynamics under physical simulation show that the system does not naturally loop back.By optimizing using gradient descent and fast projection, the dynamics can be forced to loop back, as shown in (a).The fast projection method shows better convergence in both terms of the loss functions, as shown in (b) and (c).In (d) we show our Gauss-Newton optimization with the three initial guesses: fast projection converges to the lowest loss in this case.
smooth and more natural than the gradient descent.From the good initial guess produced by fast projection, our Gauss-Newton solver converges to better local minimum Fig. 11 (d).Interestingly, if we ignore the initial state constraint in our formula by setting  = +∞, the final solution is a nice energy-preserving animation.
5.4.2Twisting bar.As mentioned in Remark 2, fast projection tends to fail when the required control force is large.Here we show an elastic bar releasing from a twisted initial state (Fig. 12).In this challenging large deformation scenario, fast projection quickly diverges, and gradient descent also fails to make progress.However, our Gauss-Newton step can still generate seamless cyclic animation using forward simulation as the initial guess.In this case, the fast projection method fails to enforce the cyclic constraint, and as such is not a good choice for the initial guess.By using the physical simulation as the initial guess, the Gauss-Newton solver generates a cyclic animation (bottom row).

LIMITATIONS AND FUTURE WORK
Self-collision.Although we integrated our method with a penalty-based contact model to include collision objects in simulations, we have not attempted to address self-collisions in our method.
Damping.Both our continuous and discrete physical systems are always Lyapunov-stable due to their cyclic nature.However, since we do not put constraints on the kinetic energy of the system, it is possible that the system arrives at a more damped solution compared to the initial, non-cyclic trajectory.A different loss function that enables more control can potentially resolve the problem.Scalability.Our Gauss-Newton solver requires us to invert the Hessian matrix (21) with size ( ×  ) 2 .While our current implementation uses a direct Cholesky decomposition for inversion, it is likely that a more efficient solver can be designed to solve a larger Hessian matrix by exploiting its structure.

CONCLUSION
We propose a constraint-free formulation for physical cyclic animations, allowing us to synthesize seamless looping animation without solving non-linear constraints.Our formulation is general and can handle a variety of physical systems.

Fig. 5 .
Fig.5.Fast projection finds a point on the constraint manifold { = 0} that is close to the initial guess u (0) = 0.The method iteratively linearizes the constraint function and project the iterand u ( ) to the zeros ( ) .
Fig. 7. Ribbons.A rotating point handle is driving three ribbons.The handle orbits in a circular motion and circles 3 times in 100 frames.The forward simulation produces an end state (Frame 99) that significantly differs from the initial state (Frame 0).
Fig.8.Bouncing Doughnut.Our method can be applied to soft body simulation with contact.Our initial guess generates an animation that almost loops back, at the cost of losing some momentum due to a large control force, making it not bounce as high.Our optimization produces a seamless loop, while preserving the momentum.

Fig. 9 .Fig. 10 .
Fig. 9.The vertical position and velocity of the center of mass of the soft body in Fig. 8 plotted against time.The forward simulation exhibits a discontinuous jump at time 0. Our initial guess (fast projection) reduces the jump, but loses the vertical momentum.After the Gauss-Newton optimization, our method restores the momentum, while maintaining a seamless physical cyclic animation.

Fig. 11 .
Fig. 11.Swinging cloth.A piece of cloth swings back and forth from a horizontal initial state.The dynamics under physical simulation show that the system does not naturally loop back.By optimizing using gradient descent and fast projection, the dynamics can be forced to loop back, as shown in (a).The fast projection method shows better convergence in both terms of the loss functions, as shown in (b) and (c).In (d) we show our Gauss-Newton optimization with the three initial guesses: fast projection converges to the lowest loss in this case.

Fig. 12 .
Fig.12.Twisting bar.An elastic bar is initialized in a twisted state.The forward simulation gradually returns to the rest state (top row).In this case, the fast projection method fails to enforce the cyclic constraint, and as such is not a good choice for the initial guess.By using the physical simulation as the initial guess, the Gauss-Newton solver generates a cyclic animation (bottom row).

Table 1 .
Performance and statistics.