How to access model jacobian from FMU or Dymola without analytical jacobian

I am trying to find a way to access the jacobian for a model in dymola either through a compiled FMU or from the exported Dymola source code.

The final objective is to use the same procedure to access the jacobian for a much more complex multibody vehicle model (205 states).

Using fmi2GetDirectionalDerivative() from the FMI Standard seemed promising so I made a simple linear vehicle model to test this.

model Vehicle "Single-track Linear bicycle vehicle model"
  extends Modelica.Blocks.Icons.Block;
  import SI = Modelica.SIunits;
  import MB = Modelica.Mechanics.MultiBody;

  // model parameters
  parameter SI.Velocity u = 10 "forward velocity";
  parameter SI.Inertia Iz = 2000 "yaw moment of inertia";
  parameter SI.Length L = 3 "wheel base";
  parameter SI.Mass Mf = 900 "front axle mass";
  parameter SI.Mass Mr = 600 "rear axle mass";
  parameter Real Cf(unit="N/rad") = 300000 "front axle cornering stiffness";
  parameter Real Cr(unit="N/rad") = 200000 "rear axle cornering stiffness";

  // calculated parameters
  final parameter SI.Mass M = Mf + Mr "mass";
  final parameter SI.Length a = Mr/Mf*L "CG position front";
  final parameter SI.Length b = L - a "CG position front";

  input SI.Angle delta "steering angle" annotation(Dialog(group="Inputs"));

public 
  SI.Velocity v "lateral velocity";
  output SI.Acceleration ay "lateral acceleration";
  SI.AngularVelocity r "yaw rate";

equation 

  ay = der(v) + u*r;
  M*(der(v) + u*r) = Cf*(delta-(v+a*r)/u) + Cr*(-(v-b*r)/u);
  Iz*der(r) = a*Cf*(delta-(v+a*r)/u) - b*Cr*(-(v-b*r)/u);

end Vehicle;

This model has:

• states - v and r
• inputs - delta
• outputs - ay

For this test,

• delta=amp*sin(2*Modelica.Constants.pi*freq*time) with

• amp = 1*Modelica.Constants.pi/180

• freq = 0.5

• Version: Dymola 2020x

• Solver: RKFIX2

• Timestep: 0.01s

• Co-simulation FMU

Since this is a linear model, the jacobian should be a constant value throughout the simulation. For this model, when I set the flag Advanced.GenerateAnalyticJacobian = true, I get the following values for the model jacobian computed from fmi2GetDirectionalDerivative() for all combinations of knowns and unknowns. In all cases, dvKnown = 1 for the function.

These values are correct based on the state space equation:

+--------------+----------+
| Derivative   | Value    |
+--------------+----------+
| der(v)/delta | 200      |
+--------------+----------+
| ay/delta     | 200      |
+--------------+----------+
| der(r)/delta | 300      |
+--------------+----------+
| der(v)/v     | -33.3333 |
+--------------+----------+
| ay/v         | -33.3333 |
+--------------+----------+
| der(r)/v     | -20      |
+--------------+----------+
| der(v)/r     | -36.6667 |
+--------------+----------+
| ay/r         | -26.6667 |
+--------------+----------+
| der(r)/r     | -70      |
+--------------+----------+

However, if I set the flag Advanced.GenerateAnalyticJacobian = false, I get completely junk values below:

+--------------+-----------+
| Derivative   | Value     |
+--------------+-----------+
| der(v)/delta | -1.57E+11 |
+--------------+-----------+
| ay/delta     | -1.57E+11 |
+--------------+-----------+
| der(r)/delta | 1.52942   |
+--------------+-----------+
| der(v)/v     | -9.12E+08 |
+--------------+-----------+
| ay/v         | -9.12E+08 |
+--------------+-----------+
| der(r)/v     | 14999.8   |
+--------------+-----------+
| der(v)/r     | 5.47E+11  |
+--------------+-----------+
| ay/r         | 5.47E+11  |
+--------------+-----------+
| der(r)/r     | -2.25E+07 |
+--------------+-----------+

I expect the value to be different from the analytical value since its numerically calculated but I don't get why its completely wrong.

I tried enabling some other flags (Advanced.AllowNumericDifferentiation, Advanced.AutomaticDifferentiation) and changing the solver to CVODE, DASSL etc. but the values remain incorrect.

Unfortunately, Dymola can't calculate an analytical jacobian for the large model so I can't use that option. All the literature I read point to fmi2GetDirectionalDerivative().

I would appreciate any inputs on how to get the model jacobian out of the FMU.

If there are other methods that can be used through Dymola, that would also work since we have a source code export license.

Cannot comment, so here is an answer that is not an answer:

The results are not complete junk: For identical values in the first table (e.g. 200 for the first two rows) you get identical values in the second table (-1.57E+11). An exception is der(v)/r and ay/r, which are identical in the second table, but maybe that is because the values are truncated.

Ask Dymola to check their implementation of fmi2GetDirectionalDerivative() in combination with Advanced.GenerateAnalyticJacobian = false.

For future reference, I was able to partially solve this issue.

I get meaningful values for the numerical jacobian (very close to the analytical values) when I export a Model Exchange FMU.

I am guessing this is some kind of bug in Dymola's Co-Simulation implementation w.r.t numerical jacobians.