Unimplemented derivatives for built-in Modelica operators

[Rüdiger Franke]

Thank you for prompt fix to #3048. FMU export works now for my models that use max! Are there more unimplemented derivatives?

translateModelFMU(Rem, version="2.0") for the model:

model Rem
  input Real x;
  input Real y;
  output Real r1, r2;
equation
  r1 = rem(x,y);
  r2 = x - div(x,y)*y;
end Rem;

results in

Error: Internal error BackendDAEOptimize.generateSymbolicJacobian failed
Error: Internal error BackendDAEOptimize.createJacobian failed

Once the reason is known, the use of wrapper functions seems to provide a workaround:

function myrem
  input Real x, y;
  output Real r;
algorithm
  r := rem(x,y);
end myrem;

function mydiv
  input Real x, y;
  output Real d;
algorithm
  d := div(x,y);
end mydiv;

model Rem
  input Real x;
  input Real y;
  output Real r1, r2;
equation
  r1 = myrem(x,y);
  r2 = x - mydiv(x,y)*y;
end Rem;

[Rüdiger Franke]

The error message does not show up anymore in r24671. Instead the dependency is just ignored in the ModelStructure of an exported FMU 2. This is of course worse than an error message.

Shouldn't Differentiate.differentiateCallExp1Arg and Differentiate.differentiateCallExpNArg at least have some else branch that raises an error -- or even better forward the differentiation to some default implementation?

Note: the above workaround generates the right ModelStructure. The calls to rem and div are wrapped into functions -- couldn't the else branch behave as if unknown calls were such user functions?

[Rüdiger Franke]

Just checked out from svn and built omc myself. The error message was back for translateModelFMU(Rem, version="2.0"). Can it be that OMEdit ignores error messages?

Anyway, I added two more rules to Differentiate.differentiateCallExpNArg, attempting to assume that the derivatives of div and rem are the same as for /. See svn diff:

Compiler/BackEnd/Differentiate.mo

@@ -1685 +1685 @@
       then
         (DAE.IFEXP(DAE.CALL(Absyn.IDENT("noEvent"),{DAE.RELATION(e1,DAE.LESS(tp),e2,-1,NONE())},DAE.callAttrBuiltinBool), res1, res2), funcs);
 
+    case ("div", {e1,e2}, DAE.CALL_ATTR(ty=tp), _, _, _, _)
+      then
+        differentiateBinary(DAE.BINARY(e1, DAE.DIV(tp), e2), inDiffwrtCref, inInputData, inDiffType, inFunctionTree);
+
+    case ("rem", {e1,e2}, DAE.CALL_ATTR(ty=tp), _, _, _, _)
+      then
+        differentiateBinary(DAE.BINARY(e1, DAE.DIV(tp), e2), inDiffwrtCref, inInputData, inDiffType, inFunctionTree);
+
   end match;
 end differentiateCallExpNArg;

This makes the FMU generation work and the generated ModelStructure is correct.

Lacking pull requests and not really knowing what I did here, can someone check if this solves the ticket?

[Martin Sjölund]

I doubt it solves anything. div and rem are not continuous-time expressions (they trigger events when they change). So the derivative of them should be 0 if anything, right?

[Adrian Pop]

Definitely the case for "rem" is wrong, it should be something like e1 - (div(e1, e2) * e2);.

[Rüdiger Franke]

Zero is right for div as it only changes at event time points. Then the existing implementation of pre appears wrong though (function differentiateCalls in file Differentiate.mo, line 1333ff):

    case (e as DAE.CALL(path=Absyn.IDENT(name = "pre")), _, _, _, _)
      then
        (e,  inFunctionTree);

(probably another ticket)

So should the implementation for event triggering mathematical functions look something like:

Compiler/BackEnd/Differentiate.mo

@@ -1573 +1573 @@
         (exp_2, funcs) = differentiateExp(exp_1, inDiffwrtCref, inInputData, inDiffType, inFuncs);
       then
        (exp_2, funcs);
+
+    // diff(ceil(exp)) = 0 because changes trigger events
+    case ("ceil", _, _, _, _, _)
+      then
+        (DAE.RCONST(0.0), inFuncs);
+
+    // diff(floor(exp)) = 0 because changes trigger events
+    case ("floor", _, _, _, _, _)
+      then
+        (DAE.RCONST(0.0), inFuncs);
+
+    // diff(integer(exp)) = 0 because changes trigger events
+    case ("integer", _, _, _, _, _)
+      then
+        (DAE.RCONST(0.0), inFuncs);
   end match;
 end differentiateCallExp1Arg;
 
@@ -1685 +1700 @@
       then
         (DAE.IFEXP(DAE.CALL(Absyn.IDENT("noEvent"),{DAE.RELATION(e1,DAE.LESS(tp),e2,-1,NONE())},DAE.callAttrBuiltinBool), res1, res2), funcs);
 
+    // diff(div(e1, e2)) = 0 because changes trigger events
+    case ("div", {e1, e2}, DAE.CALL_ATTR(ty = tp), _, _, _, _)
+      then
+        (DAE.RCONST(0.0), inFunctionTree);
+
+    // diff(mod(e1, e2)) = diff(e1 - floor(e1 / e2) * e2)
+    case ("mod", {e1, e2}, DAE.CALL_ATTR(ty = tp), _, _, _, _)
+      then
+        differentiateExp(DAE.BINARY(e1, DAE.SUB(tp), DAE.BINARY(DAE.CALL(Absyn.IDENT("floor"), {e1, e2}, inAttr), DAE.MUL(tp), e2)), inDiffwrtCref, inInputData, inDiffType, inFunctionTree);
+
+    // diff(rem(e1, e2)) = diff(e1 - div(e1, e2) * e2)
+    case ("rem", {e1, e2}, DAE.CALL_ATTR(ty = tp), _, _, _, _)
+      then
+        differentiateExp(DAE.BINARY(e1, DAE.SUB(tp), DAE.BINARY(DAE.CALL(Absyn.IDENT("div"), {e1, e2}, inAttr), DAE.MUL(tp), e2)), inDiffwrtCref, inInputData, inDiffType, inFunctionTree);
+
   end match;
 end differentiateCallExpNArg;

[Rüdiger Franke]

Or do event triggering functions probably need two implementations: one for continuous integration -- zero for div -- and another one for event updates -- division for div?

[Martin Sjölund]

der(noEvent(div())) needs an implementation separate from der(div()) in my opinion.

[Rüdiger Franke]

I'm having the use case of connected FMUs with algebraic loops in mind. This was a major motivation for introducing ModelStructure and directional derivatives in FMI 2.0.

Shouldn't this also work during event update?

For now it would be great if event triggering functions would pass FMI export and result in the correct ModelStructure.

[Rüdiger Franke]

Do the DAE.CALL's possibly have an attribute saying if in event mode or not?

Then the implementation could look like:

    case ("div", {e1, e2}, DAE.CALL_ATTR(ty = tp), _, _, _, _)
      then
        (if inEventUpdate then 
           differentiateBinary(DAE.BINARY(e1, DAE.DIV(tp), e2), inDiffwrtCref, inInputData, inDiffType, inFunctionTree)
         else
           DAE.RCONST(0.0), inFunctionTree);

[Martin Sjölund]

That would be making things too easy for implementations ;)

[Willi Braun]

Added changes in r24768 and fixed.

Replying to rfranke:

Shouldn't this also work during event update?

Yes, and it does but the events are triggered
by the functions (div, mod, e.g) and than the
originally value changes, but the derivative
is still zero also at the event, so we don't
need to take care therefore at this point.