Constrained dynamics about anything out of spring systems In - PowerPoint PPT Presentation

Penalty methods In principle, you can make just Constrained dynamics about anything out of spring systems In practice, you can make just about anything as long as it’s jello A simple example Penalty constraints A bead on a wire Why not use a spring to hold the bead on the wire? The bead can slide freely along the wire, but cannot come off it Problems: no matter how hard you pull it. weak springs won’t do the job How do we simulate the strong springs give you stiff motion of the bead when systems arbitrary forces applied to it?

First order world The real world N N In this world, f = m v In the real world, f = m a f f ˆ f N f v What is the legal velocity? What is the legal acceleration? f � f T a depends on both and N What is the legal force? v ˆ f ? the faster you’re going, the ˆ Add constraint force to cancel f faster you have to turn the illegal part of f ˆ f � Compute such that only f f � = f + ˆ f generates legal acceleration f � = f + ˆ f = − N · f f ˆ N · NN Constraint force Constraints Implicit: C ( x ) = | x | − r = 0 x All we need is calculate the constraint force θ Parametric: In order to calculate constraint force, we need to r know what legal acceleration is � cos θ � x = r sin θ

Legal acceleration Legal conditions What is the legal position? If we start with legal position and velocity C = 0 C = 0 C ( x ) = 1 2 x · x − 1 ˙ ˙ C = 0 C = 0 2 C ( x ) = 0 ¨ ¨ C = 0 C = 0 What is the legal velocity? ˙ C ( x ) = 0 ˙ C ( x ) = x · ˙ x = 0 We need only ensure the legal What is the legal acceleration? acceleration ¨ ¨ C ( x ) = 0 C ( x ) = ¨ x · x + ˙ x · ˙ x = 0 Constraint force Constraint force Use the legal condition to compute the constraint force ˙ ∂ C f = − ∂ C C ∂ x · ˆ ∂ x · f − m ∂ x ˙ x Rewrite the legal condition in a general form C ( x ) = ∂ ˙ How many variables do we have? C x + ∂ C ∂ C ¨ : constraint gradient ∂ x · ˙ ∂ x · ¨ x = 0 ∂ x x = f + ˆ C = ∂ ˙ ∂ x · f + ˆ f C x + ∂ C f Substitute ¨ x with ¨ ¨ ∂ x · ˙ = 0 m m ˙ ∂ C f = − ∂ C C Need one more condition to solve the constraint force ∂ x · ˆ ∂ x · f − m ∂ x ˙ x

Virtual work Constraint force Constraint force is passive - no energy gain or loss x | ∂ C x · ˆ ˙ f = 0 , ∀ ˙ ∂ x · ˙ x = 0 must point in the direction of ∂ C T = 1 ˆ Kinetic energy of the system: f 2 m ˙ x · ˙ x ∂ x f = λ∂ C ˆ ˆ Work done by and : ˙ f x · ˆ f T = m ˙ x · ¨ = ˙ x · f + ˙ f ∂ x x ˙ ∂ C f = − ∂ C C ˆ Substituting for in ∂ x · ˆ f ∂ x · f − m ∂ x ˙ x Make sure does no work for every legal velocity: ˆ f ∂ x · f − m ∂ ˙ x | ∂ C λ = − ∂ C C C ( x ) = ∂ C ∂ x · ˙ x x · ˆ ˙ ∂ x · ˙ ˙ f = 0 , ∀ ˙ x = 0 ∂ x · ˙ x = 0 ∂ C ∂ x · ∂ C ∂ x The Bead Example Geometric view N N f = λ∂ C f f Differentiating gives a normal vector ˆ C ∂ x ˆ v f N = ∂ C ∂ x · f − m ∂ ˙ λ = − ∂ C C x ∂ x · ˙ ∂ x ∂ x · ∂ C ∂ C ∂ x ˆ f ? When the system is at rest, ∂ x · f − m ∂ ˙ λ = − ∂ C C ∂ x · ˙ x ∂ C ∂ x · ∂ C ∂ x f = − N · f ˆ N · NN What about the case when is nonzero? ˙ x

Geometric view Feedback ∂ x · f − m ∂ ˙ λ = − ∂ C C ∂ x · ˙ x N In principle, ensuring legal acceleration can keep a ∂ C ∂ x · ∂ C ∂ x the particle exactly on the circle a x In practice, two problems cause the particle to drift x = f + ˆ ˙ f = f f · N N · ˙ x ¨ m N · NN − N · NN m − m numerical errors can accumulate when the ODE is not solved exactly ˙ N · ˙ x x = a x − ¨ N · NN constraints might not be met initially A feedback term handles both problems: C = − k s C − k d ˙ ¨ C instead of ¨ C = 0 Tinkertoys Constrained particle systems Particles: points in phase space Multiple constraints Now we know how to simulate a bead on a wire Each is a function C i ( x 1 , x 2 , . . . ) Legal state: C i ( x 1 , x 2 , . . . ) = 0 , ∀ i Apply the simple idea, we can create a constrained particle system Constraint force: linear combination of constraint gradients ∂ C i ∂ x , ∀ i

Particle system notations Constraint notations Additional notations General case ˙ C λ J J Q W ∂ C 1 ∂ C 1 ∂ C 1 ∂ ˙ ∂ ˙ ∂ ˙ q C 1 C 1 C 1 C 1 λ 1 ∂ q 1 ∂ q 2 ∂ q 3 n ∂ q 1 ∂ q 2 ∂ q 3 n λ 2 1 ∂ C 2 ∂ ˙ x 1 f 1 C 2 C 2 q : 3 n long state vector I ∂ q 1 1 ∂ q 1 m 1 x 2 f 2 I ... ... m 2 Q : 3 n long force vector ... . . . ∂ C m ∂ ˙ C m C m λ m ∂ q 3 n ∂ q 3 n : 3 n 3 n inverse mass matrix 1 W × C : m long constraint vector f n I x n m n λ : m Lagrangian multipliers J : m 3n Jacobian matrix × J : time derivative of Jacobian matrix ˙ Constraint equations Multiple force constraints To solve for force constraints, we need to solve for this C = ∂ ˙ C x + ∂ C C = ˙ ¨ ¨ J ˙ q + J ¨ q ∂ x · ˙ ∂ x · ¨ x linear system JWJ T λ = − ˙ J ˙ q − JWQ C = ∂ ˙ ∂ x · f + ˆ C x + ∂ C f C = ˙ ¨ q + JW ( Q + ˆ ¨ J ˙ Q ) = 0 ∂ x · ˙ = 0 m One nice property of is that it is symmetric JWJ T positive definite matrix ˙ ∂ C f = − ∂ C C JW ˆ Q = − ˙ ∂ x · ˆ J ˙ q − JWQ ∂ x · f − m ∂ x ˙ x Once the linear system has been solved, the vector is λ ˙ ∂ C ∂ x · λ∂ C ∂ x = − ∂ C C JWJ T λ = − ˙ multiplied by to produce the global constraint force J T ∂ x · f − m x J ˙ q − JWQ ∂ x ˙ vector ˆ Q

Constraints Parametric Constraints � cos θ Implicit: � x = r sin θ f C ( x ) = | x | − r = 0 x x Constraint is always met exactly θ θ Parametric: r r 1 degree of freedom: � cos θ � x = r sin θ Solve for ¨ θ First order world The real world In this world, f = m v In the real world, f = m a N N x = f + ˆ f T = ∂ x x = T ˙ ˙ ˙ θ x x ∂θ m T = ∂ x x = T ˙ θ = f + ˆ ˙ θ f x = ˙ T ˙ θ + T ¨ ∂θ ¨ m θ θ θ = f + ˆ f T ˙ ˆ θ = T · f f T · ˙ T ˙ θ + T · T ¨ m + T · m m ˆ θ = T · f f f = λ∂ C T · T ˙ m + T · ˆ ∂ x = λ N m θ = T · f m − T · ˙ T ˙ θ ¨ T · f T · T θ = 1 ˙ T · T m

Particle system notations Constraint notations Additional notations General case ˙ J J Q M q : 2 n long state vector q u ∂ q 1 ∂ q 1 ∂ ˙ ∂ ˙ q 1 q 1 x 1 θ 1 f 1 m 1 I ∂ u 1 ∂ u n ∂ u 1 ∂ u n ∂ q 2 ∂ ˙ q 2 x 2 θ 2 f 2 m 2 I : n long state vector ∂ u 1 ∂ u 1 u J : 2n n Jacobain matrix × . . ... . . . Q : 2 n long force vector . . . . ... ... J : time derivative of J ˙ : 2 n 2 n mass matrix M × θ n f n x n m n I ∂ q 2 n ∂ q 2 n ∂ ˙ ∂ ˙ q 2 n q 2 n ∂ u 1 ∂ u n ∂ u 1 ∂ u n Constraint equations Parametric vs. implicit θ = f + ˆ Parametric Implicit f q = M ( ˙ u ) = Q + ˆ x = ˙ T ˙ θ + T ¨ J ˙ u + J ¨ M ¨ Q ¨ m θ = T · f J T MJ ¨ u = − J T M ˙ u + J T Q JWJ T λ = − ˙ J T M ( ˙ u ) = J T Q J ˙ J ˙ q − JWQ T · ˙ T ˙ θ + T · T ¨ J ˙ u + J ¨ m J = ∂ C θ = T · f m − T · ˙ T ˙ J = ∂ q θ J T MJ ¨ u = − J T M ˙ u + J T Q where where ¨ J ˙ ∂ q ∂ u T · T T = ∂ x J = ∂ q where where Lagrangian dynamics ∂θ ∂ u

Parametric constraints Impress your friends Advantages: Fewer degrees of freedom The requirement that constraints not add or remove energy is called the Principle of Virtual Work Constraints are always met The � ’ s are called Lagrangain Multipliers Disadvantages: The derivative matrix J is called the Jacobian Matrix Hard to formulate constraints Hard to combine constraints How do we implement all this? How do we hook this up? We have a basic particle system which main job is to perform derivative evaluations We have a global matrix equation To add constraints to the basic system, we need to We want to build a model on the fly modify Each constraint function knows how to evaluate the the data structure of the basic system function itself and its various derivatives the “deriv eval” loop

Basic particle system Constrained system solver system interface solver system x 1 x 2 x n v 1 v 2 v n particles 6n particles GetDim ... f 1 f 2 f n n n time time m 1 m 2 m n x 1 x 2 x n Get/Set . . . v 1 v 2 v n State ... forces forces F 1 F 2 F p v 1 v 2 v n ... f 1 f 2 f n Deriv Eval constraints C 1 C 2 C m . . . m 1 m 2 m n A Constraint Jacobian matrix x p x q v p v q Each constraint contributes n blocks f p f q one or more blocks to the m p m q C code that evaluates Jacobian matrix ∂ ˙ ∂ C C p len C i ˙ m blocks C C ∂ q ∂ q apply_fun Sparsity: many empty blocks ˙ ˙ C C J J Modularity: let each constraint global structure compute its own blocks Constraint and particle indices ˙ ˙ x j x k ˙ λ ˆ C C J J q Q determine block location