A Symbolic Approach to the Projection Method Nam Pham Mark - PowerPoint PPT Presentation

A Symbolic Approach to the Projection Method Nam Pham Mark Giesbrecht University of Waterloo, Canada The Tenth Asian Symposium on Computer Mathematics Beijing 2012 (Pham & Giesbrecht) ASCM 2012 1 / 21

Outline Introduction 1 Constrained mechanical system The projection method Problem definition Our symbolic-numeric code-generating algorithm 2 Experimental results 3 (Pham & Giesbrecht) ASCM 2012 2 / 21

How to do a simulation of a physical mechanical system? Create a model of the system 1 Generate equations to describe the dynamic of the model 2 Solve the equations to determine the system response 3 Slider Crank Mechanism and Parallel Robot (Pham & Giesbrecht) ASCM 2012 3 / 21

Dynamics and kinematics of constrained mechanical system Kinematic constraint equations C ( x , t ) = 0 , (1) with m nonlinear algebraic equations of n generalized coordinates x 1 , · · · , x n ( m < n ). System dynamics x + C T M ¨ J λ = F , (2) where C J is the m × n Jacobian of the constraint matrix C M is an n × n symmetric generalized mass matrix λ is the ( m × 1) Lagrange multiplier Solving these DAEs for x ( t ) and λ ( t ) is computationally expensive (Pham & Giesbrecht) ASCM 2012 4 / 21

Dynamics and kinematics of constrained mechanical system Kinematic constraint equations C ( x , t ) = 0 , (1) with m nonlinear algebraic equations of n generalized coordinates x 1 , · · · , x n ( m < n ). Symbolic computation: Allow parameters z 1 , . . . , z ℓ ∈ R System dynamics x + C T M ¨ J λ = F , (2) where C J is the m × n Jacobian of the constraint matrix C M is an n × n symmetric generalized mass matrix λ is the ( m × 1) Lagrange multiplier Solving these DAEs for x ( t ) and λ ( t ) is computationally expensive Solving these systems with parameters is extremely expensive! (Pham & Giesbrecht) ASCM 2012 4 / 21

The Projection Method Blajer’s (1992) projection method: hide algebraic equations from the dynamic equations: Find a null space basis D , an n × r matrix, such that C J D = 0 or D T C T J = 0 , (3) x + C T J λ = F by D T Multiply both sides of M ¨ D T M ¨ x = D T F , (4) Now we have ODEs in x and u , which can be easily solved to determine the coordinates x , velocity u , and constraint reaction λ during simulation ˙ x = Du , (5) u = D T ( F − M ˙ D T MD ˙ Du ) , (6) λ = ( CM − 1 C T ) − 1 C ( M − 1 F − ˙ Du ) (7) (Pham & Giesbrecht) ASCM 2012 5 / 21

Numeric vs. Symbolic Modelling and Simulation Numeric • Numerical matrices are used to describe the system at a given instant in time. • Values must be given for all parameters, even if they aren’t really known. • The model must be rebuilt at every time step during simulation. (Pham & Giesbrecht) ASCM 2012 6 / 21

Numeric vs. Symbolic Modelling and Simulation Numeric Symbolic • Numerical matrices are used • All equations of motion are to describe the system at a formulated once instead of given instant in time. every step during simulation • Values must be given for all • Engineers can view the parameters, even if they governing equations in a aren’t really known. meaningful form • The model must be rebuilt at • Arbitrary substitutions for every time step during unknown quantities are not simulation. needed. (Pham & Giesbrecht) ASCM 2012 6 / 21

Numeric vs. Symbolic Modelling and Simulation Numeric Symbolic • Numerical matrices are used • All equations of motion are to describe the system at a formulated once instead of given instant in time. every step during simulation • Values must be given for all • Engineers can view the parameters, even if they governing equations in a aren’t really known. meaningful form • The model must be rebuilt at • Arbitrary substitutions for every time step during unknown quantities are not simulation. needed. Computer Algebra in Industrial Simulation MapleSim – symbolic physical modelling and simulation tool Talk tomorrow: Symbolic Computation Techniques for Advanced Mathematical Modelling by Junlin Xu (Pham & Giesbrecht) ASCM 2012 6 / 21

Our problem: Code generation for symbolic null spaces Formal definition Input : A ∈ R ( z 1 , z 2 , · · · , z ℓ ) m × n , with m ≤ n and rank r , Output : straight-line code which takes parameters α 1 , . . . , α ℓ ∈ R and evaluates a specific (consistent) basis of the null space of A : w 1 ( α 1 , . . . , α ℓ ) , w 2 ( α 1 , . . . , α ℓ ) , . . . , w n − r ( α 1 , . . . , α ℓ ) ∈ R n (Pham & Giesbrecht) ASCM 2012 7 / 21

Our problem: Code generation for symbolic null spaces Formal definition Input : A ∈ R ( z 1 , z 2 , · · · , z ℓ ) m × n , with m ≤ n and rank r , Output : straight-line code which takes parameters α 1 , . . . , α ℓ ∈ R and evaluates a specific (consistent) basis of the null space of A : w 1 ( α 1 , . . . , α ℓ ) , w 2 ( α 1 , . . . , α ℓ ) , . . . , w n − r ( α 1 , . . . , α ℓ ) ∈ R n Difficulties A is condensed with complex multivariate function Symbolic manipulation can lead to massive expression swell (Pham & Giesbrecht) ASCM 2012 7 / 21

Our problem: Code generation for symbolic null spaces Formal definition Input : A ∈ R ( z 1 , z 2 , · · · , z ℓ ) m × n , with m ≤ n and rank r , Output : straight-line code which takes parameters α 1 , . . . , α ℓ ∈ R and evaluates a specific (consistent) basis of the null space of A : w 1 ( α 1 , . . . , α ℓ ) , w 2 ( α 1 , . . . , α ℓ ) , . . . , w n − r ( α 1 , . . . , α ℓ ) ∈ R n Difficulties A is condensed with complex multivariate function Symbolic manipulation can lead to massive expression swell Previous proposed solutions Apply linear graph theory to reduce the number of equations (McPhee 2004) Fraction-free factoring to control the generation of large expression (Zhou, 2004) (Pham & Giesbrecht) ASCM 2012 7 / 21

Our problem: Code generation for symbolic null spaces Formal definition Input : A ∈ R ( z 1 , z 2 , · · · , z ℓ ) m × n , with m ≤ n and rank r , Output : straight-line code which takes parameters α 1 , . . . , α ℓ ∈ R and evaluates a specific (consistent) basis of the null space of A : w 1 ( α 1 , . . . , α ℓ ) , w 2 ( α 1 , . . . , α ℓ ) , . . . , w n − r ( α 1 , . . . , α ℓ ) ∈ R n Difficulties A is condensed with complex multivariate function Symbolic manipulation can lead to massive expression swell Advantages of our approach Very fast Partial and incremental symbolic evaluation (Pham & Giesbrecht) ASCM 2012 7 / 21

Example: Planar (2D) Slider Crank Mechanism Planar Slider Crank Mechanism with 1 degree of freedom   L 1 cos θ + L 2 sin β − s  = 0 C = L 1 sin θ − L 2 cos β − s  θ − f ( t )   − L 1 sin θ L 2 cos β − 1 C J = δ ( C ) δ ( θ, β ) = L 1 cos θ L 2 sin β 0   1 0 0 (Pham & Giesbrecht) ASCM 2012 8 / 21

Example: Spatial (3D) Slider Crank Mechanism In a slightly more complicated Spatial (3D) Slider Crank Mechanism, the second column is:   − L 2 cos ( β )   C J [ ∗ , 2 ] = − L 2 sin ( β ) cos ( α ) cos ( θ ) − L 2 sin ( β ) sin ( α ) sin ( θ )     L 2 sin ( β ) cos ( α ) sin ( θ ) − L 2 sin ( β ) sin ( α ) cos ( θ ) (Pham & Giesbrecht) ASCM 2012 9 / 21

Example: Spatial (3D) Slider Crank Mechanism 1 + x 2 , cos ( α ) = 1 − x 2 2 x 1 + x 2 where x = tan ( α Substitute sin ( α ) = 2 ) : − L 2 · 1 − x 32   1 + x 32   ( 1 − x 22 ) x 3 ( 1 − x 12 )  x 2 x 1 x 3  − 2 L 2 · ( 1 + x 22 )( 1 + x 32 )( 1 + x 12 ) − 8 L 2 · J [ ∗ ; 2 ] =   ( 1 + x 22 )( 1 + x 32 )( 1 + x 12 )     x 2 x 3 ( 1 − x 12 ) ( 1 − x 22 ) x 3 x 1   4 L 2 · ( 1 + x 22 )( 1 + x 32 )( 1 + x 12 ) − 4 L 2 · ( 1 + x 22 )( 1 + x 32 )( 1 + x 12 ) (Pham & Giesbrecht) ASCM 2012 10 / 21

Our algorithm Sketch of our approach Computing the null space using LU decomposition in a hybrid symbolic-numeric fashion Choose the ordering of row and column interchanges using 1 “indicative” numerical values Perform a symbolic LU decomposition of the “permuted” A without 2 pivoting Generate straight-line code to evaluate a null space basis at any 3 setting of the parameters (Pham & Giesbrecht) ASCM 2012 11 / 21

Algebraic static pivot selection Strategy for pivot selection Choose “random” values α 1 , . . . , α ℓ of parameters z 1 , . . . , z ℓ from 1 a finite subset S ⊆ C ; Return P , Q such that P · A ( α 1 , . . . , α ℓ ) · Q has an 2 LU -decomposition (without pivoting), using Gaussian Elimination with complete row/column pivoting. I.e., just record the row/column pivot selection. Good news: the probability of success is high (Schwarz-Zippel Lemma) Bad news: Choosing random points might be be numerically unstable... (Pham & Giesbrecht) ASCM 2012 12 / 21