Transcript 04-148ls_pres - Texas A&M University

AN EMBEDDED FUNCTION TOOL FOR MODELING
AND SIMULATING ESTIMATION PROBLEMS IN
AEROSPACE ENGINEERING
D. Todd Griffith, James D. Turner, 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
•Overview of automatic differentiation by OCEA
•Development of higher-order GLSDC algorithms
•Examples
Texas A&M University, Department of Aerospace Engineering
Introduction (1)
Estimation of dynamical systems broadly addresses the following
challenge:
“Given your best model of a dynamical system, including sensors, and
given the erroneous measured output of sensors, now what?!” (*)
Dynamics model:
Measurement model:
Typically, first order
differential correction
method utilized:
x (t )  f (t , x(t )); x0  x(t0 )
y  h( x (t ))
xˆ k 1  xˆ k  xˆ k
xˆ k  ( H TWH ) 1 H TW yk
H
h
x (t0 )
xˆ k
yk  y  h( xˆ k )
(*) Junkins First Sermon on Estimation of Dynamical Systems
Texas A&M University, Department of Aerospace Engineering
Introduction (2)
In this work, we develop higher-order (first through fourth-order)
generalizations of the Gaussian Least Squares Differential Correction
(GLSDC) algorithm.
These generalizations are implemented by using a Reversion of Series
Solution:
dx
1 d 2x
1 d 3x
1 d 4x
x  xguess    xguess 



 ...
2
3
4
ds g 2! ds g 3! ds g 4! ds g
dx
1
   G  h( xguess )
ds s 1
d 2x
dx
dx
1
2



G

G


 
ds 2 s 1
ds s 1 ds s 1
, and so on……..
Texas A&M University, Department of Aerospace Engineering
Presentation Outline
•Introduction
•Overview of automatic differentiation by OCEA
•Development of higher-order GLSDC algorithms
•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, for second-order:
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
– Automatic generation and integration of equations of
motion (AAS 04-242)
• In this paper, we consider higher-order GLSDC
algorithms………………….
Texas A&M University, Department of Aerospace Engineering
Presentation Outline
•Introduction
•Overview of automatic differentiation by OCEA
•Development of higher-order GLSDC algorithms
•Examples
Texas A&M University, Department of Aerospace Engineering
Higher-order GLSDC Algorithms (1)
First, we look at computing the required sensitivities for the Series
Reversion Solution.
Given the generic
nonlinear
measurement
model...
y  h( x (t ))
… the task is to compute
partial derivatives of this
model w.r.t. the unknown
parameters ---> Sensitivities
1st order:
 xo h  hi , j  hi , s xs , j  hi , s  sj (t , t0 )
2nd order:
2xo h  hi , jk  hi ,s xs , jk  hi ,st xt ,k xs , j
 hi , s  sjk (t , t0 , t0 )  hi , st  tk (t , t0 ) sj (t , t0 )
Need to compute first and higher-order state transition matrices……..
Texas A&M University, Department of Aerospace Engineering
Higher-order GLSDC Algorithms (2)
Notation for higher-order state transition matrices
1st order:
2nd order:
3rd order:
4th order:
ij  ij (t , to ) 
xi (t )
x j (to )
 ijk
 2 xi (t )
  ijk  t , to , to  
x j (to )xk (to )
 ijkl
 3 xi (t )
  ijkl  t , to , to , to  
x j (to )xk (to )xl (to )
 ijklm
 4 xi (t )
  ijklm  t , to , to , to , to  
x j (to )xk (to )xl (to )xm (to )
Texas A&M University, Department of Aerospace Engineering
Higher-order GLSDC Algorithms (3)
Higher-order state transition matrix differential equations
t
xi (t )  xi (to )   fi ( , x ( ))d
xi (t )  fi (t , x(t ))
to
Repeated differentiation of the integral form of the equations of motion,
followed by time differentiation, leads to first and higher-order state
transition matrix differential equation:
ij  fi ,s sj
ijk  fi ,s sjk  fi ,st tk sj
ijkl  fi ,s sjkl  fi ,su ul sjk  fi ,st tk sjl  fi ,st tkl sj  fi ,stu ul tk  sj
Fourth-order equations too long to show here, but are computed by
continuing from above………………….
Texas A&M University, Department of Aerospace Engineering
Higher-order GLSDC Algorithms (4)
Summary
The automatic differentiation capability allows the analyst to produce a
general estimation tool because the dynamical model and the measurement
model can be simply specified and differentiated. The work involved in
computing and validating partial derivatives by hand is not required.
Dynamics
model partials
Measurement
model partials
ij  fi ,s sj
ijk  fi ,s sjk  fi ,st tk sj
hi , j  hi , s  sj
hi , jk  hi , s  sjk  hi , st  tk  sj
Texas A&M University, Department of Aerospace Engineering
Presentation Outline
•Introduction and previous work
•Overview of automatic differentiation by OCEA
•Development of higher-order GLSDC algorithms
•Examples
Texas A&M University, Department of Aerospace Engineering
Ballistic Projectile Identification Problem (1)
Here, we estimate model parameters for pitch, q, and yaw, y, angle
models for an aerodynamically and inertially symmetric projectile.
Models for the angles given by:
q (t , x)  k1e t cos(1t  1 )  k2e t cos(2t   2 )
1
2
k3e3t cos(3t   3 )  k4
y (t , x )  k1e t sin(1t  1 )  k2e t sin(2t   2 )
1
2
k3e3t sin(3t   3 )  k5
14 unknown parameters include:
x   k1 , k2 , k3 , k4 , k5 , 1, 2 , 3 , 1, 2 , 3 , 1,  2 ,  3 
Texas A&M University, Department of Aerospace Engineering
Ballistic Projectile Identification Problem (2)
SUBROUTINE NONLINEAR_FX( T, EB_VAR, EB_FCTN )
USE EB_HANDLING
IMPLICIT NONE
! ARGUMENT LIST VARIABLES
REAL(DP)::T
TYPE(EB),
DIMENSION(NV), INTENT(IN ):: EB_VAR
TYPE(EB),
DIMENSION(NF), INTENT(INOUT):: EB_FCTN
! DEFINE LOCAL + EMBEDDED VARIABLES
TYPE(EB):: K1, K2, K3, K4, K5, LAM1, LAM2, LAM3, OMEG1, OMEG2, OMEG3,
DEL1, DEL2, DEL3
! ASSIGN LOCAL VARIABLES
K1=EB_VAR(1);K2=EB_VAR(2);K3=EB_VAR(3);K4=EB_VAR(4);K5=EB_VAR(5)
LAM1=EB_VAR(6);LAM2=EB_VAR(7);LAM3=EB_VAR(8)
OMEG1=EB_VAR(9);OMEG2=EB_VAR(10);OMEG3=EB_VAR(11)
DEL1=EB_VAR(12);DEL2=EB_VAR(13);DEL3=EB_VAR(14)
q
y
! COMPUTE NONLINEAR FUNCTION USING EMBEDDED ALGEBRA
EB_FCTN(1) = K1*EXP(LAM1*T)*COS(OMEG1*T+DEL1) + K2*EXP(LAM2*T)*&
COS(OMEG2*T+DEL2) + K3*EXP(LAM3*T)*COS(OMEG3*T+DEL3) + K4
EB_FCTN(2) = K1*EXP(LAM1*T)*SIN(OMEG1*T+DEL1) + K2*EXP(LAM2*T)*&
SIN(OMEG2*T+DEL2) + K3*EXP(LAM3*T)*SIN(OMEG3*T+DEL3) + K5
END SUBROUTINE NONLINEAR_FX
Texas A&M University, Department of Aerospace Engineering
Ballistic Projectile Identification Problem (3)
Cost for Ballistic Projectile Identification Problem
Iteration
1
2
3
4
5
6
7
First order
0.46453E+04
0.70828E+03
0.15229E+03
0.12652E+02
0.88894E-01
0.78899E-05
0.19504E-11
Second order
0.46453E+04
0.40676E+03
0.21104E+02
0.39243E-02
0.19141E-11
Texas A&M University, Department of Aerospace Engineering
Planar orbit example (1)
x2
Drag model:
V
r
EQM:
q
x1
 x1 (t0 )  1000 
 x (t )   2500 
 2 0  

Truth:  x1 (t0 )    200 
 x (t )   100 
 2 0  

 p   0.001 
f drag   p V V
x3
 x1  
x  
x4
 2  
x   x3     px3 V
 x   g  px V
4
 4 
 p  
0







r  
 x12  x22 

Measurements: h     


1
q  
 tan ( x2 / x1 ) 

Observed for 20 seconds at 1 second intervals
Texas A&M University, Department of Aerospace Engineering
Planar orbit example (2)
Case I: Zero drag
Convergence Study for Case I
Initial guess
0.9 xtrue , 1.0 xtrue
0.9 xtrue , 0.9 xtrue
0.9 xtrue , 0.8 xtrue
0.9 xtrue , 0.7 xtrue
0.9 xtrue , 0.6 xtrue
0.9 xtrue , 0.5 xtrue
0.9 xtrue , 0.4 xtrue
0.9 xtrue , 0.3 xtrue
0.9 xtrue , 0.2 xtrue
0.9 xtrue , 0.1 xtrue
0.9 xtrue , 0.0 xtrue
First Order
Iteration Count
5
5
5
Second Order
Iteration Count
4
4
4
5
5
4
4
6
6
4
5
6
7
5
5
7
8
5
6
Texas A&M University, Department of Aerospace Engineering
Planar orbit example (3)
Case II: With drag
Convergence Study for Case II
Initial guess
0.9 xtrue , 1.0 xtrue , p = 0.95*ptrue
0.9 xtrue , 0.9 xtrue , p = 0.95*ptrue
0.9 xtrue , 0.8 xtrue , p = 0.95*ptrue
0.9 xtrue , 0.7 xtrue , p = 0.95*ptrue
0.9 xtrue , 0.6 xtrue , p = 0.95*ptrue
0.9 xtrue , 0.5 xtrue , p = 0.95*ptrue
First Order
Iteration
Count
5
5
5
5
6
6
Second Order
Iteration
Count
4
5
6
6
7
8
Texas A&M University, Department of Aerospace Engineering
Another look at the Reversion of
Series Solution
dx
1 d 2x
1 d 3x
1 d 4x
x  xguess    xguess 



 ...
2
3
4
ds g 2! ds g 3! ds g 4! ds g
First-order
Second-order
And, so on…..
Texas A&M University, Department of Aerospace Engineering
Conclusions
•Introduced an automatic differentiation (AD) tool OCEA
•Presented higher-order Gaussian Least Squares
Differentiatal Correction (GLSDC) algorithms
•Simulated two examples using first and second-order
algorithms
•Future work includes:
–Simulation of third and fourth-order algorithms
–Extensions in propagation of uncertainty
Texas A&M University, Department of Aerospace Engineering