TestInitialConditions.py test fails on Windows

[seanmcleod]

This test fails while comparing the initial condition specified for phi in reset00.xml and the value returned from fdm['ic/phi-deg'] for the mk82_script.xml.

In script mk82_script.xml: phi should be 0.000000 but found 90.000000

The initial conditions in reset00.xml are:

<phi unit="DEG">         0.0000  </phi>
<theta unit="DEG">       90.0000  </theta>
<psi unit="DEG">         0.0000  </psi>

Placing a breakpoint in the python test these are the results when inspecting the 3 attitude angles:

fdm['ic/phi-deg']  90.0  
fdm['ic/psi-deg']  0.0  
fdm['ic/theta-deg']  90.0

I then set a breakpoint in FGInitialCondition::Load_v1() with JSBSim executing the mk82_script.xml.

After reading in the XML elements phi, theta and psi:

vOrient - {0.00000000000000000, 1.5707963267948966, 0.00000000000000000}

Which is as expected. The Euler angles from vOrient are then converted into a Quarternion.

orientation = FGQuaternion(vOrient);

orientation - 
data: {0.707, 0.000, 0.707, 0.000}
mCacheValid: false  
mEulerAngles: {0.000, -0.000, 0.000}

Now I execute the following in the debugger:

orientation.GetEulerDeg(ePhi)  90.000000000000000
orientation.GetEulerDeg(eTht)  90.000000000000000
orientation.GetEulerDeg(ePsi)  90.000000000000000

And orientation now:

orientation - 
data: {0.707, 0.000, 0.707, 0.000}
mCacheValid: true  
mEulerAngles: {1.5708, 1.5708, 1.5708}

So a couple of questions.

Why the apparent difference between the test passing on Linux (as per Travis logs) and failing on Windows?

Why the difference between the Python results for the attitude angles and the results I get in the debugger from C++ with all 3 now 90 degrees?

Is this a singularity case? But even so why the differences?

[seanmcleod]

Refreshing my basic/limited quaternions knowledge I guess I can see that the results I saw with the C++ debugger are a valid set of Euler angles (90, 90, 90) for the quaternion, but I don't see how the results from the Python test case (0, 90, 90) are.

[agodemar]

According to Euler angle theory, a value of theta = 90 deg is not allowed in the 3-2-1 convention. It's a singularity case, when the fuselage is aligned with the local vertical and nose pointing upwards.

As you see, for such a given attitude, with a given position of the wings w.r.t. Earth axes, if you assign 90 deg to theta, there are infinite possible combinations of (phi, psi). So, there are infinite triplets of Euler angles with an elevation angle theta relaxed to a value of 90 degrees.
On the contrary, this attitude is well represented by a unique (Euler-Rodriguez normalized) quaternion value.

[bcoconni]

@seanmcleod Let's keep the Python thing apart first. Do you get the same results if you execute the script directly with the JSBSim executable ? On my Linux box (fedora 27), when I run the script with the following command

$ JSBSim scripts/mk82_script.xml  --logdirectivefile=data_output/position.xml

then I get the following data in position.csv

Time	Altitude ASL (ft)	Altitude AGL (ft)	Phi (deg)	Theta (deg)	Psi (deg)
0	6594.4881989993	6594.4881989993	0	90	0

Do you get the same result from the C++ executable ?

In my experience, the results of some computations is sensitive to math optimizations made by the compiler. Do you get the same result if you disable all the optimizations from the compiler ?

[seanmcleod]

@agodemar nice diagram! I see Jon mentions briefly some of the issues with a pitch attitude of 90 degrees in the JSBSim manual in terms of the rocket launch example in terms of rocket guidance, and avoiding some of the issues etc.

@bcoconni I was wondering if some of the differences between the Windows build and the Linux build could be due to some slight differences in terms of some of the trigonometric functions like atan etc. But as you say it may also depend on compiler optimisations.

In my Python testing I was using a release version of JSBSim, since I previously had issues building a debug version of JSBSim and getting the debug version to build/link with the Python libs.

For my reports in the original post above regarding the C++ output I was seeing I was using a debug build of JSBSim running the mk82_script directly (no Python test case code involved) and reporting based on inspecting/executing code in the C++ debugger.

So I was using different compiler optimisations between the two. Also I wasn't looking at output logged, rather the variables in the debugger while executing Load_v1().

What I'll do is log the output as per your example and try both a release and debug version and report back with the results.

[bcoconni]

Some further thoughts: I just reviewed the code that converts quaternions to Euler angles and I think that the code at stake is

src/math/FGMatrix33.cpp lines 166-169

  if (data[8] == 0.0)
    mEulerAngles(1) = 0.5*M_PI;
  else
    mEulerAngles(1) = atan2(data[7], data[8]);

I guess that the initial conditions in the script mk82_script.xml corresponds to a case where both data[7] and data[8] are very close to 0.0

1. If data[8]==0.0 then the Euler angle phi is set to π/2 (this is what you get on Windows ?)
2. If the numerical round off is such that data[8] is slightly different from 0.0 then the second case of the condition applies and set phi to 0.0 because data[7]==0.0 (this is what I get on Linux ?)

[seanmcleod]

So I ran both the Windows x64 release and debug builds.

.\Release\JSBSim.exe .\scripts\mk82_script.xml --logdirectivefile=data_output\position.xml
.\Debug\JSBSim.exe .\scripts\mk82_script.xml --logdirectivefile=data_output\position.xml

The data in position.csv is identical between the release and debug builds, so no differences in terms of compiler optimisations between the two Windows builds.

Time	Altitude ASL (ft)	Altitude AGL (ft)	Phi (deg)	Theta (deg)	Psi (deg)
0	6594.4881989993	6594.4881989993	180	90	180

So what is interesting is the difference I see in the data logged to the .csv file and the data I see when I set a breakpoint in FGInitialCondition::Load_v1() and then execute orientation.GetEulerDeg(ePhi).

[bcoconni]

As @agodemar mentioned we are here in a "gimbal lock" where there is not a unique set of Euler angles to describe the position.

I am pretty sure that we are hitting a round off problem. The result of the formula I have shown above is summed up in the matrix below for a number of scenarios where data[7] and data[8] are either equal or very close to zero.

data[7]	data[8]	phi (deg)
0.0	0.0	90
0.0	1E-8	0
0.0	-1E-8	180

It illustrates that for very small variations of data[8] the angle phi has very different values. If it is confirmed that it is indeed the cause of the problem then I am unsure about the way to proceed:

1. Rewrite the Euler angles computations so that they are not so sensitive to round off (this is theoretically impossible but we might be able to alleviate the problem in most cases)
2. Force the Euler angles to special/arbitrary values for the gimbal lock so that we can predict the output for the script mk82_script
3. Modify TestInitialConditions.py to not check phi (or psi) when theta is equal to 90 deg.

[seanmcleod]

At the time when I initially came across this issue I tested a couple of values for theta in reset00.xml to see how sensitive the test was to values around 90 degrees. For example a value of 89.9 for theta will result in the test passing.

I'd vote for option 3 if it is confirmed this issue is related to the round off example you've presented.

[seanmcleod]

I ran the debugger and took a look at the values for data[7] and data[8] at the time that the output code requested the Euler angles to output them to the .csv file. And sure enough as you suspected:

data[7]	data[8]
0.00000000000000000	-2.2204460492503131e-16

Which ties up with the value of 180 degrees for phi that I see in my .csv file.

[seanmcleod]

@bcoconni, @agodemar should I go ahead and submit a commit with option 3 above implemented?

[bcoconni]

@seanmcleod Sure, go ahead.

[seanmcleod]

@bcoconni I submitted a pull request for this - #43

[bcoconni]

PR #43 has been merged (commit 470fdce)