Local, symbolic sensitivity analysis of Euler Angles

Specific Domains
1

Hey all :slight_smile:

I’m trying to wrap my head around a supposedly simple problem, but it’s not working yet.

It is known (see https://www.nxp.com/files-static/sensors/doc/app_note/AN3461.pdf) that the relationship between the measured (e.g. by an accelerometer) acceleration of earth’s gravity a=\left[a_x,a_y,a_z\right]^T relates to Euler ZYX angles (Euler angles - Wikipedia) via

\begin{bmatrix} a_x\\a_y\\a_z \end{bmatrix} = \begin{bmatrix} \sin\theta\\-\cos\theta \sin\phi\\-\cos\theta \cos\phi\end{bmatrix} (1).

\theta is the so called Pitch angle and \phi the Roll angle. They can be calculated from (1) via

\begin{eqnarray} \tan\phi&=\frac{-a_y}{-a_z}\quad (2) \\\tan\theta&=\frac{a_x\cos\phi}{-a_z} \quad (3) \\ &=\frac{a_x}{\sqrt{a_y^2+a_z^2}} \quad (4) \end{eqnarray}

What I’d like to calculate is a local sensitivity function for J_\phi(\phi,\theta), i.e. I’d like to plot a surface plot with ||{J}|| on the z-axis and \phi and \theta on the x and y axis respectively.
To my understanding, I’d need to calculate the partial derivatives of (2). The problem here is, that (2) is no function of \phi or \theta but of a_y and a_z. So there are more calculation steps required. However, (hopefully) all information are shown in (1).

Can Julia do these term conversion to eliminate all a_{x,y,z} terms? I do not know if that’s mathematically possible. But then I would not know how to calculate the sensitivity function.

I tried to use (3) as a starting point and calculated the partial derivatives, but there are still dependencies of a_z and a_x.

Can Julia do (in some magic way) a sensitivity analysis of (1) directly?

From a Monte Carlo analysis (sampling a randomly across the unit sphere’s surface and applying (2) and (3) to calculate the angles) and practical experiments, some shape like here (with different scaling on z) would be expected:
image

Disclaimer: I don’t know if that task is solvable analytically, but some paper are doing a sensitivity analysis without showing the steps (and then applying simplification like setting a_x=0 to plot it), e.g. Autocalibration of MEMS Accelerometers | IEEE Journals & Magazine | IEEE Xplore

Thank you very much!
Jan

2

Sure, you just make a grid of theta and phi values and then calculate the norm of J for each grid point and plot it

3

Thanks :slight_smile: The problem: I don’t have J(\cdot) yet and I’m wondering if Julia can calculate that via Symbolic.jl or similar tools.

4

What exactly is J? You have not defined it yet.

5

If J is the Jacobian and you don’t really need its symbolic expression, why don’t you just use ForwardDiff.jl?

1 Like
6

Sorry, it is indeed misleading in my original post. Can’t edit anymore unfortunately.

J shall be the sensitivity function for my angle \phi. Something like J_\phi(\phi,\theta)=\frac{\partial f(\phi,\theta)}{\partial\phi}. The issue is that I don’t know a formulation for f(\phi,\theta), i.e. f being a function of the angles only, and not of a_{x,y,z}. f is currently unknown to me. It should follow from (1) or (2) and (3).

My first question is if Julia can use (1) to calculate that missing f(\phi,\theta) via Symbolics.jl. I can’t do it manually, maybe it is not possible.

Second question would then be to how to create the sensitivity function J. Either from f(\phi,\theta) (if that expression can be found, then sensitivity is a straight forward partial derivative) via Symbolics.jl or ForwardDiff.jl or directly from (1) but I don’t know how to formulate that problem in Symbolics.jl or ForwardDiff.jl.

7

I think you want the magnitude of the projection of the gravitational acceleration onto an orientation vector for the accelerometer.

Also your equation (1) is missing a values on the right hand side.

Let a coordinate system be oriented with z along the gravitational direction. Then [0,0,g] is the accel vector. Let the accelerometer measure acceleration along a unit length vector [a,b,c]. Then the measured acceleration is the dot product of the two. [0,0,g] \cdot [a,b,c]

You can express [a,b,c] in terms of Euler angles, and then you still take the dot product… The sensitivity function you seek is the partial of the dot product wrt the angle of interest.

The dot product is just g*c using this coordinate system.