Constrained dynamics Simple particle system In principle, you can - PowerPoint PPT Presentation

Constrained dynamics

Simple particle system • In principle, you can make just about anything out of spring systems • In practice, you can make just about anything as long as it’s jello

Hard constraints

• Constraint force • Single implicit constraint • Multiple implicit constraint • Parametric constraint (advanced) • Implementation

A simple example A bead on a wire The bead can slide freely along the wire, but cannot come off no matter how hard you pull it. How do we simulate the motion of the bead when arbitrary forces applied to it?

Penalty constraints • Why not use a spring to hold the bead on the wire? • Problems: • weak springs won’t do the job • strong springs give you stiff systems

Constraint force • A better way would be directly computing the required constrained force that cancels only the part of applied force that acts against the constraint. applied force = sum of all external forces

First order world N f In this world, f = m v ˆ f N f What is the legal velocity? f ′ f T What is the legal force? ˆ Add constraint force to f cancel the illegal part of f f ′ = f + ˆ f f = − N · f ˆ N · NN

The real world N f In the real world, f = m a v What is the legal acceleration? a depends on both and N v ˆ f ? the faster you’re going, the faster you have to turn ˆ Compute such that only f ′ f generates legal acceleration f ′ = f + ˆ f

Constraint force • Need to compute constraint forces that cancel the illegal applied forces • Which means we need to know what legal acceleration is

• Constraint force • Single implicit constraint • Multiple implicit constraint • Parametric constraint (advanced) • Implementation

Constraints Implicit: C ( x ) = | x | − r = 0 x θ Parametric: r � � cos θ x = r sin θ

Legal velocity What is the legal position? C = 0 C ( x ) = 1 2 x · x − 1 ˙ 2 = 0 C = 0 ¨ C = 0 What is the legal velocity? ˙ C ( x ) = x · ˙ x = 0

Quiz Legal position C = 0 C ( x ) = 1 2 x · x − 1 ˙ 2 = 0 C = 0 ¨ C = 0 Legal velocity ˙ C ( x ) = x · ˙ x = 0 What is the legal acceleration? ¨ C ( x ) = ¨ x · x + ˙ x · ˙ x = 0

Legal conditions If we start with legal position C = 0 and velocity ˙ C = 0 ¨ C ( x ) = 0 C = 0 ˙ C ( x ) = 0 We need only to ensure the legal acceleration ¨ C ( x ) = 0

Constraint force Use the legal condition to compute the constraint force Rewrite the legal condition in a general form C ( x ) = ∂ ˙ x + ∂ C C ∂ C ¨ : constraint gradient ∂ x · ˙ ∂ x · ¨ x = 0 ∂ x x = f + ˆ f x with Substitute ¨ ¨ m C = ∂ ˙ ∂ x · f + ˆ C x + ∂ C f ¨ ∂ x · ˙ = 0 m ˙ ∂ C f = − ∂ C C ∂ x · ˆ ∂ x · f − m x ∂ x ˙

Constraint force ˙ ∂ C f = − ∂ C C ∂ x · ˆ ∂ x · f − m x ∂ x ˙ How many variables do we have? Need one more condition to solve the constraint force

Virtual work Constraint force is passive - no energy gain or loss T = 1 Kinetic energy of the system: 2 m ˙ x · ˙ x ˆ Virtual work done by and : ˙ x · ˆ f f T = ˙ x · m ¨ = ˙ x · f + ˙ f x ˆ Make sure does no work for every legal velocity: f x | ∂ C C ( x ) = ∂ C x · ˆ ˙ ∂ x · ˙ ˙ f = 0 , ∀ ˙ x = 0 ∂ x · ˙ x = 0

Constraint force x | ∂ C x · ˆ ∂ x · ˙ ˙ f = 0 , ∀ ˙ x = 0 must point in the direction of ∂ C ˆ f ∂ x f = λ∂ C ˆ ∂ x ˙ ∂ C f = − ∂ C C ˆ ∂ x · ˆ Substituting for in f ∂ x · f − m x ∂ x ˙ ∂ x · f − m ∂ ˙ λ = − ∂ C C x ∂ x · ˙ ∂ C ∂ x · ∂ C ∂ x

Two conditions • Legal acceleration ˙ ∂ C f = − ∂ C C ∂ x · ˆ ∂ x · f − m x ∂ x ˙ • Principle of virtual work f = λ∂ C ˆ ∂ x

The Bead Example N f C ( x ) = 1 2 x · x − 1 2 = 0 v f = λ∂ C ˆ ˆ f ? ∂ x ∂ x · f − m ∂ ˙ λ = − ∂ C C x ∂ x · ˙ ∂ C ∂ x · ∂ C ∂ x

Quiz x = (0 , − 2) C ( x ) = 1 2 x · x − 1 2 r 2 x = ( − 4 , 0) ˙ m = 10 kg g = (0 , − 9 . 8) f = λ∂ C r = 2 ˆ ∂ x ∂ x · f − m ∂ ˙ λ = − ∂ C C x ∂ x · ˙ ∂ C ∂ x · ∂ C ∂ x What is the constraint force at the current state?

Feedback • In principle, ensuring legal acceleration can keep the particle exactly on the circle • In practice, two problems cause the particle to drift • numerical errors can accumulate when the ODE is not solved exactly • constraints might not be met initially

Feedback • A feedback term handles both problems: C = − k s C − k d ˙ ¨ C instead of ¨ C = 0

Quiz What is the new formula for lambda after adding the feedback term? C = − k s C − k d ˙ ¨ C ∂ x · f − m ∂ ˙ λ = − ∂ C C x ∂ x · ˙ ∂ C ∂ x · ∂ C ∂ x

• Constraint force • Single implicit constraint • Multiple implicit constraint • Parametric constraint (advanced) • Implementation

Tinkertoys • Now we know how to simulate a bead on a wire • Apply the simple idea, we can create a constrained particle system

Constrained particle system • Particles: each particle stores its current position and velocity • Forces: each force affects the acceleration of certain particles • Constraints: Each is a function C i ( x 1 , x 2 , . . . ) Legal state: C i ( x 1 , x 2 , . . . ) = 0 , ∀ i Constraint force on each particle: linear combination of constraint gradients ∂ C i ∂ x , ∀ i

Constraint gradients Normal of C 1 : ∂ C 1 ∂ x Normal of C 2 : ∂ C 2 ∂ x Legal states: the intersection of two planes Normal of the legal states: ∂ C 1 ∂ C 2 λ 1 ∂ x + λ 2 ∂ x

Implicit constraint • Each constraint is represented by an implicit function • What does the normal of a hypersurface mean? The direction where the particle is not allowed to move • What does the intersection of the hypersurface represent? Legal state • The constraint force lies in the space spanned by the normals of constraints

Particle system notations General 3D case Q W q 1 f 1 x 1 I q : 3 n long position vector 1 m 1 x 2 f 2 I m 2 ... Q : 3 n long force vector . . . : 3 n 3 n inverse mass matrix W × 1 I f n x n m n

Constraint notations Additional notations ˙ λ J J ∂ ˙ ∂ ˙ ∂ ˙ ∂ C 1 ∂ C 1 ∂ C 1 C 1 C 1 C 1 λ 1 ∂ q 1 ∂ q 2 ∂ q 3 n ∂ q 1 ∂ q 2 ∂ q 3 n λ 2 ∂ ˙ ∂ C 2 C 2 ∂ q 1 ∂ q 1 ... ... ∂ ˙ ∂ C m C m λ m ∂ q 3 n ∂ q 3 n λ : m Lagrangian multipliers J : m 3n Jacobian matrix × J : time derivative of Jacobian matrix ˙

Constraint equations C = ∂ ˙ C x + ∂ C C = ˙ ¨ ¨ J ˙ q + J ¨ q ∂ x · ˙ ∂ x · ¨ x C = ∂ ˙ ∂ x · f + ˆ C x + ∂ C f C = ˙ ¨ q + JW ( Q + ˆ ¨ J ˙ Q ) = 0 ∂ x · ˙ = 0 m ˙ ∂ C f = − ∂ C C JW ˆ Q = − ˙ ∂ x · ˆ J ˙ q − JWQ ∂ x · f − m x ∂ x ˙ ˙ ∂ C ∂ x · λ∂ C ∂ x = − ∂ C C JWJ T λ = − ˙ J ˙ q − JWQ ∂ x · f − m x ∂ x ˙

Multiple force constraints To solve for force constraints, we need to solve for this linear system JWJ T λ = − ˙ J ˙ q − JWQ What is the dimension of ? JWJ T Once the linear system has been solved, the vector is λ multiplied by to produce the global constraint force J T vector ˆ Q

Impress your friends • The requirement that constraints not add or remove energy is called the Principle of Virtual Work • The λ ’ s are called Lagrangain Multipliers • The derivative matrix J is called the Jacobian Matrix

Quiz C ( x 1 , x 2 ) = 1 2( x 1 − x 2 ) · ( x 1 − x 2 ) − 1 2 d 2 What’s the Jacobian matrix of this constraint?

Quiz C ( x 1 ) = 1 2 x 1 · x 1 − 1 2 r 2 C ( x 1 , x 2 ) = 1 2( x 1 − x 2 ) · ( x 1 − x 2 ) − 1 2 d 2 What’s the Jacobian matrix for this system?

• Constraint force • Single implicit constraint • Multiple implicit constraint • Parametric constraint (advanced) • Implementation

Constraints Implicit: C ( x ) = | x | − r = 0 x θ Parametric: r � � cos θ x = r sin θ

Parametric constraints � � cos θ x = r sin θ f x Constraint is always met exactly θ r 1 degree of freedom: ¨ Solve for θ

Equation of motion N T = ∂ x x = T ˙ ˙ θ x ∂θ θ = f + ˆ f x = ˙ T ˙ θ + T ¨ ¨ m θ ˆ θ = T · f f T · ˙ T ˙ θ + T · T ¨ m + T · m T ˙ m − T · ˙ θ = T · f θ ¨ T · T

Particle system notations General 2D case Q M q : 2 n long position vector q u f 1 m 1 I θ 1 x 1 f 2 x 2 θ 2 m 2 I : n long parameter vector u ... . . . . . . Q : 2 n long force vector . . . : 2 n by 2 n mass matrix M θ n f n m n I x n

Constraint notations Additional notations ˙ J J ∂ q 1 ∂ q 1 ∂ ˙ ∂ ˙ q 1 q 1 ∂ u 1 ∂ u n ∂ u 1 ∂ u n ∂ q 2 ∂ ˙ q 2 ∂ u 1 ∂ u 1 J : 2n by n Jacobian matrix ... ... J : time derivative of J ˙ ∂ q 2 n ∂ q 2 n ∂ ˙ ∂ ˙ q 2 n q 2 n ∂ u 1 ∂ u n ∂ u 1 ∂ u n