Derivatives in initialization equations

[baggepinnen]

MTK requires that all derivatives that are referred to in initialization equations and the varmap are present also in the regular equations. In modelica, one can specify der(y) = 0 in initialization equations despite der(y) not appearing in any regular equation, which is nice because it's common to want to initialize at zero derivative of outputs of block components.

Here's a link to a simple modelica model that demonstrates this

https://playground.modelica.university/?model=ICBwYXJhbWV0ZXIgUmVhbCBlPTAuOCAiQ29lZmZpY2llbnQgb2YgcmVzdGl0dXRpb24iOwoKICBjb25zdGFudCBSZWFsIGVwcz0xZS0zICJTbWFsbCBoZWlnaHQiOwogIEJvb2xlYW4gZG9uZSAiRmxhZyB3aGVuIHRvIHR1cm4gb2ZmIGdyYXZpdHkiOwogIFJlYWwgaCAiSGVpZ2h0IjsKICBSZWFsIHYgIlZlbG9jaXR5IjsKICBSZWFsIHk7CmluaXRpYWwgZXF1YXRpb24KICBoID0gNy4wICJJbml0aWFsIGhlaWdodCI7CiAgZGVyKHkpID0gMC4wICJJbml0aWFsIGJ1dCBleHByZXNzZWQgaW4gdGVybXMgb2YgZGVyIjsKICAvL3YgPSAwLjAgIkluaXRpYWwgYnV0IGV4cHJlc3NlZCBpbiB0ZXJtcyBvZiBkZXIiOwogIGRvbmUgPSBmYWxzZTsKZXF1YXRpb24KICB2ID0gZGVyKGgpOwogIHkgPSBoOwogIGRlcih2KSA9IGlmIGRvbmUgdGhlbiAwIGVsc2UgLTkuODE7CiAgd2hlbiB7aDwwLGg8LWVwc30gdGhlbgogICAgZG9uZSA9IGg8LWVwczsKICAgIHJlaW5pdCh2LCAtZSooaWYgaDwtZXBzIHRoZW4gMCBlbHNlIHByZSh2KSkpOwogIGVuZCB3aGVuOwo%253D&report=VGhlIGNsYXNzaWMgTW9kZWxpY2EgImJvdW5jaW5nIGJhbGwiLiAgUHJlc3MgdGhlICJQbGF5IiBidXR0b24gdG8gYW5pbWF0ZS4KCkNoYW5nZSB0aGUgY29lZmZpY2llbnQgb2YgcmVzaXR1dGlvbiB0byBzZWUgaG93IGl0IGNoYW5nZXMgdGhlIGJvdW5jaW5nIG9mIHRoZSBiYWxsLgoKPHN2ZyBoZWlnaHQ9IjQwMCIgd2lkdGg9IjEwMCUiIHZpZXdib3g9IjAgMCA0MDAgNDAwIj4KICA8Y2lyY2xlIGN4PSJ7eyAxMDArbm93LnRpbWUqMjUgfX0iIGN5PSJ7eyAxMDArNyoyNSAtMjUqbm93LmggfX0iIHI9IjQwIiBzdHJva2U9ImJsYWNrIiBzdHJva2Utd2lkdGg9IjMiIGZpbGw9InJlZCIgLz4KICA8cmVjdCB4PSIwIiB5PSJ7ezE0MCs3KjI1fX0iIHdpZHRoPSI0MDAiIGhlaWdodD0iMTAiIHN0cm9rZT0iYmxhY2siLz4KPC9zdmc%252BCgo8Y2hhcnQgc2lnbmFscz0iaCI%252BPC9jaGFydD4%253D

The model is reproduced below in case the link stops working

  parameter Real e=0.8 "Coefficient of restitution";

  constant Real eps=1e-3 "Small height";
  Boolean done "Flag when to turn off gravity";
  Real h "Height";
  Real v "Velocity";
  Real y;
initial equation
  h = 7.0 "Initial height";
  der(y) = 0.0 "Initial vel but expressed in terms of der of output";
  //v = 0.0 "Initial vel";
  done = false;
equation
  v = der(h);
  y = h;
  der(v) = if done then 0 else -9.81;
  when {h<0,h<-eps} then
    done = h<-eps;
    reinit(v, -e*(if h<-eps then 0 else pre(v)));
  end when;

[AayushSabharwal]

structural_simplify should be able to recognize that D(y) is present in the system and handle it appropriately. To avoid having to specify constraints that are binding as initialization_eqs are, we can also search defaults. The user can remove the initial condition without having to re-simplify the system. For a nicer API, structural_simplify could take a higher_order_derivatives keyword which takes an array of higher order derivatives to compute.

[AayushSabharwal]

I'm not fully sure about this, but maybe we can handle this outside of structural simplification too? If y is an observed variable y ~ f(..) then D(y) ~ D(f(..)) and we should be able to use expand_derivatives to reduce this. Recursively applying this procedure to the RHS and combining with dummy_sub for known derivatives should just work?

[matthew-kapp]

Is initialization_eqs available as a macro yet?

[ChrisRackauckas]

If it's not working right now (I thought it was?), then you'd just recurse the derivative map https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/nonlinear/initializesystem.jl#L44-L47 until the derivatives are all gone. Basically structural simplify gives a matching of derivative expressions to equations, and then once you have that you can just recurse that to get the basic form of any expression, and that can then be used to generate the non-derivative form of the initialization equation.