Issue with MultiBody models in the NF due to missing evaluation during flattening

[Francesco Casella]

After @adrpo's commits for overconstrained connectors, the MultiBody examples all flatten and produce a balanced system of equations. Howevver, there are still problems that prevents them to run.

Please check ​Modelica.Mechanics.MultiBody.Examples.Constraints.UniversalConstraint. The following error is reported:

[/var/lib/hudson/slave/workspace/OpenModelica_TEST_LIBS/OpenModelica/OMCompiler/Compiler/BackEnd/Differentiate.mo:402:3-402:155:writable] Error: Derivative of expression "joint.revolute_b.R_rel.T[1,2] = (matrix({1.0, 0.0, 0.0}) * transpose(matrix({1.0, 0.0, 0.0})))[1, 2] - (matrix({1.0, 0.0, 0.0}) * transpose(matrix({1.0, 0.0, 0.0})))[1, 2] * cos(joint.phi_b)" w.r.t. "time" is non-existent.
[/var/lib/hudson/slave/workspace/OpenModelica_TEST_LIBS/OpenModelica/OMCompiler/build/lib/omlibrary/Modelica trunk/Mechanics/MultiBody/Joints.mo:306:7-306:60:writable] Error: Internal error 
Differentiate.differentiateEquationTime failed for joint.revolute_b.R_rel.T[1,2] = (matrix({1.0, 0.0, 0.0}) * transpose(matrix({1.0, 0.0, 0.0})))[1, 2] - (matrix({1.0, 0.0, 0.0}) * transpose(matrix({1.0, 0.0, 0.0})))[1, 2] * cos(joint.phi_b)
Error: Internal error - IndexReduction.pantelidesIndexReduction1 failed! Use -d=bltdump to get more information.

For some reason, the expression

(matrix({1.0, 0.0, 0.0}) * transpose(matrix({1.0, 0.0, 0.0})))[1, 2] 
- (matrix({1.0, 0.0, 0.0}) * transpose(matrix({1.0, 0.0, 0.0})))[1, 2] * cos(joint.phi_b)

is not evaluated appropriately during flattening, causing a failure of Pantelides's algoritm later on.

[Per Östlund]

I'm not sure why a different error is reported in the error log compared to the one in this ticket, but the actual error is currently:

Error: Derivative of expression "joint.revolute_b.R_rel.T[1,1] = (matrix(joint.revolute_b.e) * transpose(matrix(joint.revolute_b.e)))[1, 1] + (1.0 - (matrix(joint.revolute_b.e) * transpose(matrix(joint.revolute_b.e)))[1, 1]) * cos(joint.phi_b) - (skew(joint.revolute_b.e) * sin(joint.phi_b))[1, 1]" w.r.t. "time" is non-existent.

The issue seems to be in the Modelica.Mechanics.MultiBody.Frames.TransformationMatrices.planarRotation function, which is defined as:

function planarRotation
  input Real e[3];
  input Modelica.SIunits.Angle angle;
  output TransformationMatrices.Orientation T;
algorithm
  T := outerProduct(e,e) + (identity(3) - outerProduct(e,e))*Math.cos(angle) -
       skew(e)*Math.sin(angle);
  annotation(Inline = true);
end planarRotation;

The NF doesn't mess around much with algorithms, and flattens it as:

  T := matrix(e) * transpose(matrix(e)) + 
       ({{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}} -
       matrix(e) * transpose(matrix(e))) * cos(angle) - skew(e) * sin(angle);

The OF on the other hand produces:

  T := {{ e[1] ^ 2.0 + (1.0 - e[1] ^ 2.0) * cos(angle),
         (e[1] - e[1] * cos(angle)) * e[2] + e[3] * sin(angle),
         (e[1] - e[1] * cos(angle)) * e[3] - e[2] * sin(angle)},
        {(e[2] - e[2] * cos(angle)) * e[1] - e[3] * sin(angle),
          e[2] ^ 2.0 + (1.0 - e[2] ^ 2.0) * cos(angle),
         (e[2] - e[2] * cos(angle)) * e[3] + e[1] * sin(angle)},
        {(e[3] - e[3] * cos(angle)) * e[1] + e[2] * sin(angle),
         (e[3] - e[3] * cos(angle)) * e[2] - e[1] * sin(angle),
          e[3] ^ 2.0 + (1.0 - e[3] ^ 2.0) * cos(angle)}};

The issue occurs when the function is inlined by the backend, which produces expressions that it can't handle for some reason. "Simplifying" algorithms like the OF does is a slippery slope, and not always possible anyway (what if e[3] instead was e[:] for example?). So this seems like something that would be best handled by the backend, either by simplifying the inlined expressions or fixing the Differentiate module so it can handle this.

[Vitalij Ruge]

I think the Backend don't know e.g. skew, because the OF evaluate skew.

model test
function foo
  input Real[3] e;
  output Real[:, :] y = skew(e);
end foo;
  Real [:,:] y = foo({1,1,1});
equation

end test;

where Instantiation return

function test.foo
  input Real[3] e;
  output Real[3, 3] y = {{0.0, -e[3], e[2]}, {e[3], 0.0, -e[1]}, {-e[2], e[1], 0.0}};
end test.foo;

class test
  Real y[1,1];
  Real y[1,2];
  Real y[1,3];
  Real y[2,1];
  Real y[2,2];
  Real y[2,3];
  Real y[3,1];
  Real y[3,2];
  Real y[3,3];
equation
  y = {{0.0, -1.0, 1.0}, {1.0, 0.0, -1.0}, {-1.0, 1.0, 0.0}};
end test;

So, I would rule out the error in inline module.
In case e[:] I get a crash in the runtime for OF.

[Per Östlund]

Replying to vitalij:

So, I would rule out the error in inline module.

The inline module does what it's supposed to do in this case, but the inlined expression can't be handled by e.g. the differentiation even though it's a perfectly valid expression. So maybe some simplification needs to be applied after the inlining.

However, I now realised that skew does in fact always take a vector of three elements and returns a 3x3 matrix, so that's something the NF could handle by just adding the appropriate definition of skew to ModelicaBuiltin.

[Per Östlund]

Fixed in 9688d8a and 74c9f25.