{"id":112,"date":"2011-02-01T20:00:54","date_gmt":"2011-02-01T20:00:54","guid":{"rendered":"https:\/\/blogs.mathworks.com\/seth\/2011\/02\/01\/wrong-answers-and-gimbal-lock\/"},"modified":"2011-02-01T20:53:07","modified_gmt":"2011-02-01T20:53:07","slug":"wrong-answers-and-gimbal-lock","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/simulink\/2011\/02\/01\/wrong-answers-and-gimbal-lock\/","title":{"rendered":"Wrong Answers and Gimbal Lock"},"content":{"rendered":"<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/seth\/guy_rouleau_small.png\" alt=\"Blogger, Guy Rouleau\" style=\"float: left; margin-right: 1em;\">By <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/authors\/31651\">Guy Rouleau<\/a><br><br>\r\n\r\nI received a request last week where the code generated using Real-time Workshop was giving different answers from the original simulation. This type of requests is the one that scares me the most, because the goal of Real-Time Workshop is to generate a code equivalent to a Simulink model.<\/p>\r\n\r\n<p>Let's see how I approached this challenge!<\/p>\r\n\r\n<p><strong>1. Reproducing the issue<\/strong><\/p>\r\n\r\n<p>The first thing I do when I have a \u201cdifferent results\u201d case is to confirm what the user reports. I built the model using the GRT target and compared the logged data. They are different!<\/p>\r\n\r\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/seth\/2010Q4\/InitialResults.png\" alt=\"Comparaison of results and generated code\" \/><\/p>\r\n\r\n<p><strong>2. Understanding the system<\/strong><\/p>\r\n\r\n<p>The model had an architecture similar to the <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2010b\/toolbox\/aeroblks\/f4-50264.html\">NASA HL-20 demo<\/a> shipping with the <a href=\"https:\/\/www.mathworks.com\/products\/aeroblks\/\">Aerospace Blockset<\/a>.<\/p>\r\n\r\n<p>My mental representation of this model looks like:<\/p>\r\n\r\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/seth\/2010Q4\/ModelSchema.png\" alt=\"Model Schematic\" \/><\/p>\r\n\r\n<p>This architecture is used in many domains where the motion of a body in space must be simulated, for example robotics and aerospace. <\/p>\r\n\r\n<p><strong>3. Understanding the results<\/strong><\/p>\r\n\r\n<p>Let\u2019s take a closer look at the different data and zoom on where they diverge:<\/p>\r\n\r\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/seth\/2010Q4\/Data.png\" alt=\"Comparaison of results and generated code, zoomed\" \/><\/p>\r\n\r\n<p>This signal represents an orientation oscillating between +90 and -90 degrees and we see that there is a clear switch around -90 degrees.<\/p>\r\n\r\n<p><strong>4. Guessing something<\/strong><\/p>\r\n\r\n<p>This is when Seth passed by my office. Based on the indices we had, he suggested that the model might be very sensitive to numerical perturbations. It is well known that floating point round off can cause small differences on the order of <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2010b\/techdoc\/ref\/eps.html\">eps<\/a> between generated code and simulation. But when a system is very sensitive and unstable, a difference even as small as <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2010b\/techdoc\/ref\/eps.html\">eps<\/a> can grow into completely different results.<\/p>\r\n\r\n<p><strong>5. Testing the hypothesis<\/strong><\/p>\r\n\r\n<p>If the model is sensitive to numerical perturbations, modifying the initial states by a random value on the order of <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2010b\/techdoc\/ref\/eps.html\">eps<\/a> will affect the results significantly. I tested it, and it changes the simulation results significantly. What I notice is that every time I run the simulation using a new random value, the \u201cswitch\u201d around 90 degrees happens at different times.<\/p>\r\n\r\n\r\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/seth\/2010Q4\/PerturbationResults.png\" alt=\"Comparaison of results with numerical perturbation\" \/><\/p>\r\n\r\n<p><strong>6. Finding a solution<\/strong><\/p>\r\n\r\n<p>Now that we know that the model is numerically sensitive, it would be great to fix it!<\/p>\r\n\r\n<p>Having generated six degrees of freedom trajectories for robotic manipulators in the past, I remembered learning somewhere that the <a href=\"http:\/\/en.wikipedia.org\/wiki\/Euler_angles\">Euler angles<\/a> representation of orientations has a singularity, known as the <a href=\"http:\/\/en.wikipedia.org\/wiki\/Gimbal_lock\">gimbal lock<\/a>.<\/p>\r\n\r\n<p>In the <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2010b\/toolbox\/aeroblks\/6dofeulerangles.html\">6Dof Euler Angles<\/a> block, the angles are obtained through integration of the Euler Angles derivative, which is computed using this equation:<\/p> \r\n\r\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/seth\/2010Q4\/EulerAnglesEquation.gif\" alt=\"Relationship to determine the Euler rate vector\" \/><\/p>\r\n\r\n<p>where [<i>p q r<\/i>]<sup>T<\/sup> is the body-fixed angular velocity vector and <img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/seth\/2010Q4\/phi_theta_psi.gif\"\/> is the the rate of change of the Euler angles. In this equation, when <em>theta<\/em> is close to 90 degrees the matrix becomes badly scaled and causes the numerical sensitivity described above.<\/p>\r\n\r\n<p> An alternative to integrate the body motion is to compute the rate of change of <a href=\"http:\/\/en.wikipedia.org\/wiki\/Quaternions_and_spatial_rotation\">rotation quaternions<\/a> instead. :<\/p>\r\n\r\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/seth\/2010Q4\/QuaternionEquation.gif\" alt=\"Integration using quaternions\" \/><\/p>\r\n\r\n<p>Quaternions have a lot of advantages: They are numerically more stable than rotation matrices and avoid the gimbal lock problem of Euler angles.<\/p>\r\n\r\n<p><strong>7. Validate your solution<\/strong><\/p>\r\n\r\n<p>In the Simulink model I replaced the Euler Angles implementation:<\/p>\r\n\r\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/seth\/2010Q4\/EulerImplementation.png\" alt=\"Motion Integration using Euler Angles\" \/><\/p>\r\n\r\n<p>by the quaternion representation:<\/p>\r\n\r\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/seth\/2010Q4\/QuaternionImplementation.png\" alt=\"Motion Integration using quaternion\" \/><\/p>\r\n\r\n<p>I tested the solution by comparing simulation with numerical perturbation and generated code. Now simulation results and generated code results are matching as expected!<\/p>\r\n\r\n<p><strong>Now it\u2019s your turn!<\/strong><\/p>\r\n\r\n<p>What is your process of investigation when Simulink does not give the results you expect? Leave a <a href=\"https:\/\/blogs.mathworks.com\/seth\/?p=112&amp;#comment\">comment here<\/a><\/p>\r\n\r\n<p><\/p>","protected":false},"excerpt":{"rendered":"<p>By Guy Rouleau\r\n\r\nI received a request last week where the code generated using Real-time Workshop was giving different answers from the original simulation. This type of requests is the one that... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/simulink\/2011\/02\/01\/wrong-answers-and-gimbal-lock\/\">read more >><\/a><\/p>","protected":false},"author":41,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[121,43,76,33],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/simulink\/wp-json\/wp\/v2\/posts\/112"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/simulink\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/simulink\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/simulink\/wp-json\/wp\/v2\/users\/41"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/simulink\/wp-json\/wp\/v2\/comments?post=112"}],"version-history":[{"count":0,"href":"https:\/\/blogs.mathworks.com\/simulink\/wp-json\/wp\/v2\/posts\/112\/revisions"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/simulink\/wp-json\/wp\/v2\/media?parent=112"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/simulink\/wp-json\/wp\/v2\/categories?post=112"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/simulink\/wp-json\/wp\/v2\/tags?post=112"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}