Transcript AUTOMATIC GENERATION AND INTEGRATION OF EQUATIONS

AUTOMATIC GENERATION AND INTEGRATION OF
EQUATIONS OF MOTION BY OPERATOR OVER
LOADING TECHNIQUES
(x ,y )
3

3

1

2

d  L  L
 Q  CT λ
 
dt  q  q

N
D. Todd Griffith, Andrew J. Sinclair, James D. Turner,
John E. Hurtado, and John L. Junkins
Texas A&M University
Department of Aerospace Engineering
College Station, TX 77840
Texas A&M University, Department of Aerospace Engineering
Presentation Outline
•Introduction and previous work
•Overview of automatic differentiation by OCEA
•Equation of motion formulation using automatic
differentiation
•A new algorithm for numerical integration
•Examples
Texas A&M University, Department of Aerospace Engineering
Introduction and previous work
• Work on multibody dynamics dates back to around 1960’s
• Multibody dynamics is a mature field; many codes and
many books have been written in the area
– Authors include: Shabana, Scheihlen, Garcia de Jalon
• In general, many methods exist for equation of motion
generation. Typically, Lagrange’s equations, Kane’s
equations, or Newton/Euler Methods employed
• The primary question is which method for generation of
equations of motion is most suitable for problem
complexity, generality desired, computational issues, and
computational resources to name a few.
Texas A&M University, Department of Aerospace Engineering
Equation of motion formulation using
automatic differentiation
• Lagrange’s equations:
Lagrangian: L  T V
d  L  L
T
 QC λ
 
dt  q  q
subject to
Cq  b
T, V: kinetic and
potential energy
q: generalized
coordinates
Q: generalized force
C: constraint matrix
l: Lagrange multipliers
Of course, many choices for equation of motion formulation exist;
however, the utility of automatic differentiation is immediately seen for
implementing Lagrange’s Equations………………...
Texas A&M University, Department of Aerospace Engineering
Presentation Outline
•Introduction and previous work
•Overview of automatic differentiation by OCEA
•Equation of motion formulation using automatic
differentiation
•A new algorithm for numerical integration
•Examples
Texas A&M University, Department of Aerospace Engineering
Overview of automatic differentiation by OCEA(1)
• OCEA (Object Oriented Coordinate Embedding Method)
– Extension for FORTRAN90 (F90) written by James D. Turner
– Automatic Differentiation (AD) equation manipulation package in which
new data types are created in order to define independent and dependent
variables (functions)
• OCEA description
– Automatic differentiation enabled by coding rules for differentiation Chain rule of calculus
– Can compute first through fourth-order partial derivatives of scalar,
vector, matrix, and higher order tensors
– Standard mathematical library functions (e.g. sin, cos, exp) are overloaded
– Scalars are replaced by differential n-tuple
f :  f
f
 2 f 
Texas A&M University, Department of Aerospace Engineering
Overview of automatic differentiation by OCEA(2)
• OCEA description (cont’d)
– Intrinsic operators such as ( +, - , * , / , = ) are also overloaded to enable,
for example, addition and multiplication of OCEA variables:
a  b :  a  b a  b  2 a   2b 
a * b :  a * b  i  a * b   j  i  a * b  
– Partial derivative computation is hidden to the user - takes place in the
background without user intervention.
– Access to partial derivatives (e.g. Jacobian, Hessian, etc.) made simple by
overloading of the “ = “ sign.
Texas A&M University, Department of Aerospace Engineering
Overview of automatic differentiation by OCEA(3)
• The need for computation of partial derivatives is found in
numerous applications
• Previous applications include
– Root solving and Optimization
– Second and higher-order GLSDC algorithms (AAS 04148)
• In this paper, we consider automatic generation of
equations of motion for linked mechanical systems
Texas A&M University, Department of Aerospace Engineering
Presentation Outline
•Introduction and previous work
•Overview of automatic differentiation by OCEA
•Equation of motion formulation using automatic
differentiation
•A new algorithm for numerical integration
•Examples
Texas A&M University, Department of Aerospace Engineering
Approach (1): Direct approach(1)
By differentiating the Lagrangian
d  L  L
T


Q

C
λ
 
dt  q  q
d  T 
 2T
 2T
qj 
qj


dt  qi  qi q j
qi q j
 2T
M  mij 
qi q j
 2T
M  mij 
qi q j
d  T  L
T


Q

C
λ


dt  q  q
L
 Automatic Differentiation (AD)
q
Q  Specified
C
AD or specified
l
many methods
Mass matrix and its time derivative computed
by second-order differentiation………..
Texas A&M University, Department of Aerospace Engineering
Approach (1): Direct approach(2)
By differentiating the Lagrangian
Can also compute constraint matrix, C,
automatically for holonomic type constraint.
 (q )  0
 (q ) 




q
= Cq 
C 
q
t
t
q
Now forming
equations:
L
Mq + Mq 
 Q  CT λ
q
 2T
M  mij 
qi q j
 2T
M  mij 
qi q j
Accelerations computed
after generating or
prescribing all terms here.
L
-1 
T 
q = M   Mq 
+ Q  C λ  Now, can proceed with
q

 numerical integration…….
Texas A&M University, Department of Aerospace Engineering
Approach (2): A modified form of Lagrange’s
Equations
d  L  L
T


Q

C
λ
 
dt  q  q
In this approach, we can bypass
forming the kinetic energy function.
Must still form potential energy (if you
want to), but most importantly need
an efficient method for computing the
mass matrix and partial derivatives of
mass matrix…………..
V (q)
M ( q ) q  G (q , q ) 
 Q  CT λ
q
G (q, q)  qT H (1) q ... qT H ( n ) q 
H
(i )
h
(i )
jk
1  mij mik m jk
 



2  qk
q j
qi

V
T 
q = M  G 
+ Q  C λ
q


-1
Texas A&M University, Department of Aerospace Engineering



Approach (2): Identifying the mass matrix (1)
KE for system of rigid bodies
1 n
T (v , ω)    mi viT vi  I iiT i 
2 i 1
Introduce the
following velocity
transformation
matrices:
ri
vi  Ai (q)q 
q
q
 i
ωi  Bi (q)q 
q
q
For the case of planar
rigid bodies,
transformation matrices
specialize to partials of
position of mass center
and partials of angular
orientations
Additionally, Bi is a constant vector of 1’s and 0’s for a
minimal coordinate planar system
Texas A&M University, Department of Aerospace Engineering
Approach (2): Identifying the mass matrix (2)
1 n
T (v , ω)    mi viT vi  I iiT i 
2 i 1

Introduce
transformations
1 n
T (q, q)   mi qiT AiT Ai qi  I i qiT BiT Bi qi
2 i 1

1 n T
  qi  mi AiT Ai + I i BiT Bi  qi
2 i 1

1 n T
  qi M (q)qi
2 i 1



Need not specify T,
can get M by
forming mass center
positions vectors…...
n
Mass matrix and
mass matrix partials
M (q )   mi AiT Ai + I i BiT Bi
i 1
 AiT
 BiT
M (q) n
T Ai 
T Bi 
  mi 
Ai + Ai

I
B
+
B
 i

i
i
q

q

q

q

q
i 1




Texas A&M University, Department of Aerospace Engineering
Presentation Outline
•Introduction and previous work
•Overview of automatic differentiation by OCEA
•Equation of motion formulation using automatic
differentiation
•A new algorithm for numerical integration
•Examples
Texas A&M University, Department of Aerospace Engineering
Numerical Integration: Solving
xk 1  xk 
h
 k0  2k1  2k2  k3 
6
k0  f ( xk , tk )
x  f ( x (t ), t )
x (t0 )  x0
Fourth-order Runge-Kutta
algorithm
Utilizes approximate derivatives
through fourth order
k1  f ( xk  k0 2 , tk  h 2 )
k2  f ( xk  k1 2 , tk  h 2 )
k3  f ( xk  k2 , tk  h)
So called Taylor integration scheme:
x (t  t )  x (t )  f ( x (t ), t ) t 

1
1
f ( x (t ), t ) t 2  f ( x (t ), t ) t 3
2!
3!
1
1
f ( x (t ), t ) t 4  f ( x (t ), t ) t 5  .....
4!
5!
Can compute through
fourth-order time
derivatives exactly.
Potentially fifth-order
method.
Texas A&M University, Department of Aerospace Engineering
Presentation Outline
•Introduction and previous work
•Overview of automatic differentiation by OCEA
•Equation of motion formulation using automatic
differentiation
•A new algorithm for numerical integration
•Examples
Texas A&M University, Department of Aerospace Engineering
Spring Pendulum by Direct method (1)
T  12 m(r 2  r 2 2 )
r
V  Vspring  Vgrav
θ
L  T V
 12 k (r  r0 )2  mg (r0  r cos )

L
T 
q = M   Mq 
+ Q  C λ
q


-1
 2T
M  mij 
qi q j
 2T
M  mij 
qi q j
Texas A&M University, Department of Aerospace Engineering
Spring Pendulum by Direct method (2)
SUBROUTINE SPRING_PEND_EQNS( PASS, TIME, X0, DXDT, FLAG )
USE EB_HANDLING
IMPLICIT NONE
****************************************
!.....LOCAL VARIABLES
TYPE(EB)::L, T, V
! LAGRANGIAN, KINETIC, POTENTIAL
REAL(DP):: M, K
! MASS AND STIFFNESS VALUES
REAL(DP), DIMENSION(NV):: JAC_L
REAL(DP), DIMENSION(NV,NV):: HES_L
****************************************
T = 0.5D0*M*(X0(3)**2 + X0(1)**2*X0(4)**2) ! DEFINE KE
V = 0.5D0*K*(X0(1)-R0)**2 + M*GRAV*(R0-X0(1)*COS(X0(2))) ! DEFINE PE
L = T – V ! DEFINE LAGRANGIAN FUNCTION
JAC_L = L
JAC_L_Q = JAC_L(1:NV/2)
! dL/(dq,dqdot)
! dL/dq
HES_L = L ! EXTRACT SECOND ORDER PARTIALS OF LAGRANGIAN
MASS
= HES_L(NV/2+1:NV,NV/2+1:NV) ! COMPUTE MASS MATRIX
MASSDOT = HES_L(NV/2+1:NV,1:NV/2)
!COMPUTE MDOT
****************************************
DXDT(1)%E = X0(3)%E
! RDOT
DXDT(2)%E = X0(4)%E
! THETADOT
DXDT(3)%E = QDOTDOT(1)
! RDOTDOT
DXDT(4)%E = QDOTDOT(2)
! THETADOTDOT
END
SUBROUTINE
SPRING_PEND_EQNS
Texas
A&M University,
Department
of Aerospace Engineering
Open-chain Topology of Rigid Bodies (1)
Mass center location can be recursively
formulated and utilized to compute
velocity transformation matrices and
ultimately the mass matrix and its
partial derivatives.
1
2
N
Generation of Equations of motion for
open-chain topologies can be
automated to the point of simply
specifying the number of bodies and
the initial conditions, while for closed
chain-topologies we additionally need
to prescribe the constraint equations.
Texas A&M University, Department of Aerospace Engineering
Open-chain Topology of Rigid Bodies (2)
10 body model for laboratory
deployment of inflated aerospace
structure.
Model easily
updated to include
joint torques (e.g.
damping)
Texas A&M University, Department of Aerospace Engineering
Closed-chain Topology Example
5 link manipulator example.

(x3,y3)


Two holonomic constraints specified.
Driving torque on link 1.


Payload Motion
Texas A&M University, Department of Aerospace Engineering
Conclusion
•Introduced an automatic differentiation (AD) tool: OCEA
•Implemented Lagrange’s equations using AD
–Mass matrix and it’s time derivatives computed by
direct differentiation of T and by recursive formulation
•Suggested a new algorithm for numerical integration
•Examples included simulations of a generalized code for
planar n-body systems with open and closed-chain
topologies
•Future work includes computational studies, extension to
3D topologies, and efforts to validate solutions for flexible
body systems
Texas A&M University, Department of Aerospace Engineering