Improper use of local frame rotation rate #553
Comments
|
For the sake of completeness there is a second occurrence of jsbsim/src/initialization/FGInitialCondition.cpp Lines 1414 to 1419 in d9df845 |
|
The bug can be fixed by forcing --- a/src/initialization/FGInitialCondition.cpp
+++ b/src/initialization/FGInitialCondition.cpp
@@ -1191,10 +1191,7 @@ bool FGInitialCondition::Load_v1(Element* document)
// Refer to Stevens and Lewis, 1.5-14a, pg. 49.
// This is the rotation rate of the "Local" frame, expressed in the local frame.
const FGMatrix33& Tl2b = orientation.GetT();
- double radInv = 1.0 / position.GetRadius();
- FGColumnVector3 vOmegaLocal = {radInv*vUVW_NED(eEast),
- -radInv*vUVW_NED(eNorth),
- -radInv*vUVW_NED(eEast)*tan(position.GetLatitude())};
+ FGColumnVector3 vOmegaLocal = { 0., 0., 0.,};
vPQR_body = Tl2b * vOmegaLocal;
@@ -1413,10 +1410,7 @@ bool FGInitialCondition::Load_v2(Element* document)
// Refer to Stevens and Lewis, 1.5-14a, pg. 49.
// This is the rotation rate of the "Local" frame, expressed in the local frame.
- double radInv = 1.0 / position.GetRadius();
- FGColumnVector3 vOmegaLocal = { radInv*vUVW_NED(eEast),
- -radInv*vUVW_NED(eNorth),
- -radInv*vUVW_NED(eEast)*tan(position.GetLatitude())};
+ FGColumnVector3 vOmegaLocal = { 0., 0., 0. };
if (attrate_el) {@seanmcleod, @jonsberndt, @agodemar I'm waiting for your feedback before removing |
|
The formula |
In trying to get up to speed I've been reading my copy of Stevens & Lewis which happens to be the 3rd edition, I'm guessing you're quoting from the 2nd or 1st edition? What I noticed is that in the 3rd edition the angular velocity isn't split into 3 terms. So in the 3rd edition's equation 1.7-16a And the
|
Yes I'm quoting the 2nd edition because it's the only edition that I could find on the internet as I don't have an access to my own personal copy (3rd edition as well) at the moment.
That's my understanding as well.
Correct. |
|
This is only a problem at a boundary case, right? Just checking. The EOM has been verified against a few NASA simulations (several years ago), so I know that generally speaking the EOM is solid. |
Sort of.
This problem is not involving the EOM integration but the initialization of The sad thing is that |
|
So in summary, looking at the current initialization code. So if we're dealing with a V1 initialization file then angular rates for the aircraft can't be specified by the user. void FGPropagate::SetInitialState(const FGInitialCondition* FGIC)
{
...
// Set the angular velocities of the body frame relative to the ECEF frame,
// expressed in the body frame.
VState.vPQR = FGIC->GetPQRRadpsIC();
VState.vPQRi = VState.vPQR + Ti2b * in.vOmegaPlanet;
}
/** Gets the initial body rotation rate
@return Initial body rotation rate in radians/second */
const FGColumnVector3 GetPQRRadpsIC(void) const { return vPQR_body; }
bool FGInitialCondition::Load_v1(Element* document)
{
...
vPQR_body = Tl2b * vOmegaLocal;
}
bool FGInitialCondition::Load_v2(Element* document)
{
...
if (attrate_el) {
string frame = attrate_el->GetAttributeValue("frame");
frame = to_lower(frame);
FGColumnVector3 vAttRate = attrate_el->FindElementTripletConvertTo("RAD/SEC");
if (frame == "eci") {
FGMatrix33 Ti2l = position.GetTec2l() * Ti2ec;
vPQR_body = Tl2b * Ti2l * (vAttRate - vOmegaEarth);
} else if (frame == "ecef") {
vPQR_body = Tl2b * position.GetTec2l() * vAttRate;
} else if (frame == "local") {
vPQR_body = Tl2b * (vAttRate + vOmegaLocal);
} else if (frame == "body") {
vPQR_body = vAttRate;
} else if (!frame.empty()) { // misspelling of frame
cerr << endl << fgred << " Attitude rate frame type: \"" << frame
<< "\" is not supported!" << reset << endl << endl;
result = false;
} else if (frame.empty()) {
vPQR_body = Tl2b * vOmegaLocal;
}
} else { // Body frame attitude rate assumed 0 relative to local.
vPQR_body = Tl2b * vOmegaLocal;
}
}
So given the current code in JSBSim listed above in terms of initialization I'm not so sure that
Snippets from original JSBSim Reference Manual PDF.
So taking the |
|
@bcoconni @jonsberndt any further thoughts? |
|
Oh yes , many thoughts but I need to put my ideas in writing. You made a good point but as far as I'm concerned that's not the end of the story 🙂 Following your comment I'd like to change the title of this issue to Improper use of local frame rotation rate. Hopefully I'll be able to contribute something later this week. |
|
First of all, I'd like to stress that the current issue is not the result of a formula that we are trying to use outside of its domain of validity. The formula 1.5-14a from Stevens & Lewis 2nd edition is an exact formula that is giving a correct result for all latitudes from -90 degrees to +90 degrees. However... Analogy with the gimbal lock problem for Euler anglesWhat we are facing here is an issue similar to the gimbal lock for Euler angles. The conversion from matrix transformations to the Euler angles is using exact formulas, still those are giving This has been worked around in JSBSim by forcing the yaw angle to zero (which is an arbitrary choice) when the gimbal is, well, locked: jsbsim/src/math/FGMatrix33.cpp Lines 180 to 187 in d3536e9 As a consequence the Euler angles rates are undefined for the gimbal lock situation. I also would like to continue this analogy to address @jonsberndt's comment:
The issue with the local frame rotation rate is just as much a boundary case than the gimbal lock is a boundary case for the Euler angles. This does not allow us to dismiss the problem since we don't want JSBSim to crash abruptly or give silly results when the aircraft nose is pointing to the sky, do we ? Even so if the likeliness of flying straight to the sky at an exact pitch of +90 degrees is tenuous. Thus I see no reason why the current issue (undefined local frame rotational rates at the Poles) could be dismissed because it is very unlikely that someone would fly very close to the Poles. We are aiming at a robust FDM software that behaves well even in corner cases, aren't we ? Initializing
|
bcoconni
changed the title
I can't say for sure what @jonsberndt's thinking is, but my guess is he was thinking that it could potentially be solved along the lines of the gimbal lock, i.e. boundary check for a latitude of +-90 degrees and do something in this case to avoid tan blowing up etc., as opposed to simply ignoring the crash. |
The problem is that you can "do something" about undefined angles (either Euler or North/East), for instance by forcing their value to an arbitrary value. By undefined value it is meant that an infinity of values are a solution to the problem. So you pick one value among the infinite many possible solutions and bam! the problem vanishes. However you can't do anything about their rate because determining a function derivative at an undefined value is undecidable. That's the very reason why their derivatives diverge to infinity. Even Euler angles rates can't be computed at the gimbal lock so JSBSim is just keeping the previous rates unchanged i.e. using their previous finite value whatever that was: jsbsim/src/models/FGAuxiliary.cpp Lines 142 to 145 in db4d967 However this strategy somehow works because JSBSim does not do anything from the Euler rates. They are just computed for users to do whatever they see fit from these parameters. But the situation is much different for What I don't understand is why I have built a small example that shows that the current code forces a vehicle to spin in the absence of forces and moments (since we don't need to make any assumption about forces and moments to obtain that result). I'd be curious to know why that is the expected and correct behaviour ? |
|
Ok, I'm just now noticing this thread, and just waking up. I'll try and take a look at this in a couple of hours.
|
|
I'm definitely no expert in terms of the equations of motion and their implementation with respect to all the relevant reference frames, so please double-check my comments 😉 First off, are we not missing So we have PropertyManager->Tie("velocities/p-rad_sec", this, eP, (PMF)&FGPropagate::GetPQR);
PropertyManager->Tie("velocities/q-rad_sec", this, eQ, (PMF)&FGPropagate::GetPQR);
PropertyManager->Tie("velocities/r-rad_sec", this, eR, (PMF)&FGPropagate::GetPQR);
PropertyManager->Tie("velocities/pi-rad_sec", this, eP, (PMF)&FGPropagate::GetPQRi);
PropertyManager->Tie("velocities/qi-rad_sec", this, eQ, (PMF)&FGPropagate::GetPQRi);
PropertyManager->Tie("velocities/ri-rad_sec", this, eR, (PMF)&FGPropagate::GetPQRi);So if an aircraft is running down a runway due east with some velocity, then The point with Glancing at the JSBSim source it looks like the actual integration and propagation takes place in the ECI frame. So |
I'm not advocating for it be kept at all costs, I'm still getting up to speed on the issues and options 😉
I agree JSBSim should behave well even in corner cases. Although we could make this corner case a smaller corner than 'fly very close to the Poles'. It's only an issue with initialization as you mentioned, so if the aircraft was initialized away from the Pole it could still fly over/close to the Pole. Plus it's only an issue if the initialization is at or close to the Poles and some initial translational velocity is specified. In other words the case of an aircraft starting it's take-off at the Pole with an initial 0 velocity would be fine if we modified the following code to check if double radInv = 1.0 / position.GetRadius();
FGColumnVector3 vOmegaLocal = { radInv*vUVW_NED(eEast),
-radInv*vUVW_NED(eNorth),
-radInv*vUVW_NED(eEast)*tan(position.GetLatitude())}; |
There is a strong assumption behind this statement: you are assuming that the runway follows the Earth curvature. Indeed in that case, Now if we assume the same situation (a runway due East with a non zero velocity) but this time we assume that the runway has a slope then your assumption breaks down since
Of course the picture above is grossly exaggerated but the point is that you cannot assume that the aircraft is following the Earth curvature when running on a runway. Now let's assume the aircraft is initialized with some eastward velocity on the runway shown above then |
Sorry, I have overreacted there but I'd really like the discussion to focus on the reason why we are using At the moment, I am not convinced at all that initializing
In fact, it does not have to be a corner case at all since in my opinion it is the result of an incorrect initialization.
It's also an issue when a satellite is initialized with a polar orbit. In that case, the easiest way to initialize the satellite is to spawn it over the Poles with a
I am not sure I am following your reasoning here, since if we only use the formula when |
This example made me realize that using Say the satellite is initialized very close to the Pole (not exactly at the Pole because East would be meaningless) with an Eastward velocity then Now if the same satellite is initialized on the same polar orbit but just above the equator then its Eastward velocity will be zero meaning that So the initial value of |
No problem. But I was making the point/asking whether the example of an aircraft with an initial translational velocity doesn't imply that
But I can see now, based on your most recent replies that there is an underlying assumption that these So I'm leaning towards your point that
The point I was trying to make was that currently the code throws an exception from But if you as you're suggesting we get rid of |
Thanks. Let's wait for @jonsberndt opinion.
Below is a diff to ======================================================================
FAIL: testIndependenceOfInitialLocation (__main__.TestInitialConditions)
----------------------------------------------------------------------
Traceback (most recent call last):
File "/home/bertrand/OpenSource/JSBSim/GitHub/JSBSim/tests/TestInitialConditions.py", line 413, in testIndependenceOfInitialLocation
self.assertAlmostEqual(fdm['ic/r-rad_sec'], r)
AssertionError: 0.0 != -0.6315477934082796 within 7 places (0.6315477934082796 difference)
----------------------------------------------------------------------Basically, the test spawns the satellite very close to the North Pole (at a latitude of 89.9 degrees i.e. 35000 ft away on the Earth surface from the Pole) with an eastward velocity. The initial values of the As the failure message above shows, the diff --git a/tests/TestInitialConditions.py b/tests/TestInitialConditions.py
index 12710cf9..2baf1053 100644
--- a/tests/TestInitialConditions.py
+++ b/tests/TestInitialConditions.py
@@ -23,7 +23,9 @@
import os, math, shutil
import xml.etree.ElementTree as et
import pandas as pd
-from JSBSim_utils import append_xml, ExecuteUntil, JSBSimTestCase, RunTest, CreateFDM
+from JSBSim_utils import append_xml, ExecuteUntil, JSBSimTestCase, RunTest, CopyAircraftDef
# Values copied from FGJSBBase.cpp and FGXMLElement.cpp
convtoft = {'FT': 1.0, 'M': 3.2808399, 'IN': 1.0/12.0}
@@ -45,7 +47,7 @@ class TestInitialConditions(JSBSimTestCase):
IC_file = append_xml(use_tag.attrib['initialize'])
IC_tree = et.parse(os.path.join(path_to_jsbsim_aircrafts, IC_file))
- return (tree, IC_tree)
+ return tree, IC_tree, IC_file
def test_initial_conditions_v1(self):
prop_output_to_CSV = ['velocities/vc-kts']
@@ -107,7 +109,7 @@ class TestInitialConditions(JSBSimTestCase):
for s in self.script_list(('ZLT-NT-moored-1.xml',
'737_cruise_steady_turn_simplex.xml')):
- (tree, IC_tree) = self.getElementTrees(s)
+ tree, IC_tree, _ = self.getElementTrees(s)
IC_root = IC_tree.getroot()
# Only testing version 1.0 of init files
@@ -257,7 +259,7 @@ class TestInitialConditions(JSBSimTestCase):
def test_geod_position_from_init_file_v2(self):
for s in self.script_list(('ZLT-NT-moored-1.xml',
'737_cruise_steady_turn_simplex.xml')):
- (tree, IC_tree) = self.getElementTrees(s)
+ tree, IC_tree, _ = self.getElementTrees(s)
IC_root = IC_tree.getroot()
# Only testing version 2.0 of init files
@@ -355,10 +357,59 @@ class TestInitialConditions(JSBSimTestCase):
# Regression test for the bug reported by @fernandafenelon and fixed by
# Sean McLeod in #545.
def testClimbRateSetting(self):
- fdm = CreateFDM(self.sandbox)
- fdm.load_script(self.sandbox.path_to_jsbsim_file('scripts',
- 'c1722.xml'))
+ fdm = self.create_fdm()
+ self.load_script('c1722.xml')
fdm['ic/gamma-deg'] = 4
self.assertAlmostEqual(fdm['ic/gamma-deg'], 4)
+ # Regression test for issue #553
+ #
+ # This tests verifies that the initial value for the vehicle rotational rates
+ # P, Q, R do not depend on the position at which the vehicle is initialized.
+ # First reference : close to the North Pole heading East
+ # Second reference : above the Equator heading North
+ def testIndependenceOfInitialLocation(self):
+ script_path = self.sandbox.path_to_jsbsim_file('scripts/ball.xml')
+ tree, aircraft_name, _ = CopyAircraftDef(script_path, self.sandbox)
+ tree.write(self.sandbox('aircraft', aircraft_name, aircraft_name+'.xml'))
+ # Alter the initial conditions XML file to force the initial latitude
+ # to 90 degrees.
+ _, IC_tree, IC_file = self.getElementTrees(script_path)
+ IC_root = IC_tree.getroot()
+ lat_tag = IC_root.find('latitude')
+ psi_tag = IC_root.find('psi')
+ alt_tag = IC_root.find('altitude')
+ psi_tag.text = '90.0' # Heading East
+ lat_tag.text = '89.9' # Above the North Pole
+ h0 = float(alt_tag.text)
+ IC_tree.write(os.path.join('aircraft', aircraft_name, IC_file))
+
+ fdm = self.create_fdm()
+ fdm.set_aircraft_path('aircraft')
+ self.load_script('ball.xml')
+ fdm.run_ic()
+
+ p = fdm['ic/p-rad_sec']
+ q = fdm['ic/q-rad_sec']
+ r = fdm['ic/r-rad_sec']
+
+ self.delete_fdm()
+
+ # Since the equatorial radius is 70159 ft larger than the polar radius
+ # we need to decrease the altitude by the same amount in order to
+ # initialize the vehicle at the same radius.
+ alt_tag.text = str(h0-70159)
+ psi_tag.text = '0.0' # Heading North
+ lat_tag.text = '0.0' # Above equator
+ IC_tree.write(os.path.join('aircraft', aircraft_name, IC_file))
+
+ fdm = self.create_fdm()
+ fdm.set_aircraft_path('aircraft')
+ self.load_script('ball.xml')
+ fdm.run_ic()
+
+ self.assertAlmostEqual(fdm['ic/p-rad_sec'], p)
+ self.assertAlmostEqual(fdm['ic/q-rad_sec'], q)
+ self.assertAlmostEqual(fdm['ic/r-rad_sec'], r)
+ |
|
So, This is an interesting discussion. I haven't jumped into the code to look closely at this - I have some obligations today, but here are my immediate thoughts:
It sounds like you are recommending to remove vOmegaLocal. It also looks like you've done a lot of work to verify why it should be removed and that the proposed fix works. I don't remember exactly why things were done the way they were back when that code was initially written. So, I'd say, go ahead with your recommended fix. If anything in the documentation needs to be clarified or improved (or something placed on the wiki at least) that would be good to do - maybe capture some of what you wrote above and place it in the wiki. A decent portion of the JSBSim user base are engineering students, and it's good to have some kind of traceability to textbooks in common use. That's my two cents. :-) |
Sorry there again for having overreacted to your statement about the boundary case. I was afraid that this argument was used to dismiss the problem altogether. That was certainly not the case, sorry for having misinterpreted your meaning.
That would be great ! As you suggested the NESC test cases have not yet made their way to the JSBSim regression test suite. Even if your files are not perfectly tidied up that'd be a great addition to JSBSim.
That's indeed what we are trying to do with the test cases other than NESC's. The more tests, the better.
According to my investigations, all the tests succeed once
As @seanmcleod pointed out in one of his comment above, Stevens & Lewis 3rd edition does not even mention
In the short term the commits will reference this issue so that the link is maintained (at least as long as we are using GitHub). But indeed the discussion should be compiled in the wiki or in JSBSim docs.
|
Improper use of the local frame rotation rate leads to the divergence vehicle rotation rates.
Improper use of the local frame rotation rate leads to the divergence vehicle rotation rates.
Improper use of the local frame rotation rate leads to the divergence vehicle rotation rates.
Improper use of the local frame rotation rate leads to the divergence vehicle rotation rates.
|
As @seanmcleod rightly pointed out, we still have the issue of tan(90) blowing up when the frame jsbsim/src/initialization/FGInitialCondition.cpp Lines 1417 to 1425 in 68c4136 Given that in that particular case (
Any suggestion ? |
It seems to me that 3 & 4 are out. Maybe I missed this from above, but can the aircraft be initialized to some rates at the poles where the frame is not the local frame? Maybe the option is to issue a warning and gracefully exit when the (hopefully rare) condition of trying to do this occurs? Or, is there special handling that can be done just at the poles? It's early for me and I have not had my coffee yet ... |
Current frame options for the setting of initial rates. bool FGInitialCondition::Load_v2(Element* document)
{
...
...
// Initialize vehicle body rates
// Allowable frames
// - ECI (Earth Centered Inertial)
// - ECEF (Earth Centered, Earth Fixed)
// - Body
Element* attrate_el = document->FindElement("attitude_rate");
if (attrate_el) {
string frame = attrate_el->GetAttributeValue("frame");
frame = to_lower(frame);
const FGMatrix33& Tl2b = orientation.GetT();
FGColumnVector3 vAttRate = attrate_el->FindElementTripletConvertTo("RAD/SEC");
if (frame == "eci") {
FGMatrix33 Ti2l = position.GetTec2l() * Ti2ec;
vPQR_body = Tl2b * Ti2l * (vAttRate - vOmegaEarth);
} else if (frame == "ecef") {
vPQR_body = Tl2b * position.GetTec2l() * vAttRate;
} else if (frame == "local") {
// Refer to Stevens and Lewis, 1.5-14a, pg. 49.
// This is the rotation rate of the "Local" frame, expressed in the local frame.
double radInv = 1.0 / position.GetRadius();
FGColumnVector3 vOmegaLocal = {radInv*vUVW_NED(eEast),
-radInv*vUVW_NED(eNorth),
-radInv*vUVW_NED(eEast)*tan(position.GetLatitude())};
vPQR_body = Tl2b * (vAttRate + vOmegaLocal);
} else if (frame == "body") {
vPQR_body = vAttRate;
} else if (!frame.empty()) { // misspelling of frame
cerr << endl << fgred << " Attitude rate frame type: \"" << frame
<< "\" is not supported!" << reset << endl << endl;
result = false;
} else if (frame.empty()) {
vPQR_body.InitMatrix();
}
} else { // Body frame attitude rate assumed 0 relative to local.
vPQR_body.InitMatrix();
}
|
Agreed. I had some further thoughts about option 1. We should not focus on the component values of |
I'm submitting a ...
Describe the issue
JSBSim computes the rotation rate of the local frame in
src/initialization/FGInitialCondition.cpp:jsbsim/src/initialization/FGInitialCondition.cpp
Lines 1191 to 1199 in d9df845
Problem: the formula from Stevens & Lewis uses the tangent of the latitude which obviously triggers math issues when the latitude is initialized to 90 degrees.
What is the current behavior?
When the initial latitude is modified to 90 degrees in
aircraft/ball/reset00.xmlthen the tests
CheckScriptsandTestAccelerometerfail:What is the expected behavior?
The tests should run successfully. There is no reason JSBSim should fail when the latitude initialized to the legal value of 90 degrees.
What is the motivation / use case for changing the behavior?
A review of the version control history shows the following facts:
vOmegaLocalas the rotation rate of the local frame which was added toVState.vPQRjsbsim/src/models/FGPropagate.cpp
Lines 221 to 231 in a3819b9
vOmegaLocalfromsrc/models/FGPropagate.Cpptosrc/initialization/FGInitialConditions.cppHowever
VState.vPQRis measuring the angular velocity of the vehicle relative to the ECEF frame:jsbsim/src/models/FGPropagate.h
Lines 111 to 114 in d9df845
So I am challenging the relevance of adding
vOmegaLocal(rotation rate of the local frame with respect to the ECEF) tovPQR(angular velocity of the body frame with respect to the ECEF). Said otherwise, I think the vectorvOmegaLocalshould be removed altogether from JSBSim.Other information
My point is that the formula 1.5-14a from Stevens & Lewis is relevant if the inertial angular velocity is split in 3 terms:
w_b/i = w_b/v + w_v/e + w_e/i(the variablevOmegaLocalbeingw_v/eandin.OmegaPlanetisw_e/iin JSBSim parlance):However, in JSBSim the inertial angular velocity is split in 2 terms:
w_b/i = w_b/e + w_e/iwhereVState.vPWRiisw_b/i,VState.vPQRisw_b/eandin.OmegaPlanetisw_e/i. Said otherwise, there is no need for the termw_v/ei.e.vOmegaLocal.The text was updated successfully, but these errors were encountered: