Jacobians in daeMode for Dynawo

[casella]

Today I had a discussion with @marcochiaramello, @rosiereflo, and @joyelfeghali about using C code generated with --daeMode --generateSymbolicJacobian.

Consider the following test case:

model DAE
  Real u "known input variable";
  Real v1, v2, v3, v4, v5 "algebraic variables";
  Real x(start = 0, fixed = true) "state variable";
equation
  u = time;
  v1 + v2 = u;
  v1 - v2 = -u;
  v3 = v1 + v2;
  der(x) = -x + v3;
  v4 + v5 = x;
  v4 - v5 = -x;
  v6 = v4*v5;
annotation(
__OpenModelica_commandLineOptions = "--daeMode --generateSymbolicJacobian");
end DAE;

The DAE system can be expressed as

F(x, der(x), v, t) = 0   (1)

where x are the differentiated variables, v the algebraic variables, and t is the time variable. We can define the vector of the unknown of the original DAE as

z = [ x
      v ]

so the DAE can be written as

F(z, der(z), t) = 0   (2)

Assuming the DAE (2) is index-1, it can be solved numerically by IDA, which needs a function to compute the residual (2) and its Jacobians dF/dz and dF/dder(z). However, by default OMC doesn't throw the raw, original DAEs to the IDA solver, mainly for two reasons:

1. the system may be index-2 or index-3, which IDA cannot handle
2. it is advisable to keep the number of equations handled by the iterative IDA solver to the minimum amount, resorting to explicit solutions computed by the generated C code as much as possible

What OMC does is basically to re-use the causalization and index reduction algorithms that it already implements to reduce (possibly high-index) DAEs into ODEs. This allows to expose to the IDA numerical solver only the differential equations and the implicit algebraic equations that must be solved to compute the derivatives, which are matched to the state variables (after index reduction) and to the unknowns of those implicit equations. Define

z = [ x
      w ]

where x are the state variables (after index reduction) and w are the algebraic variables involved in implicit systems that are necessary to compute the derivatives.

Consider now the above-described DAE model, which is index 1, so no index reduction is required. The BLT of the system has a first 2x2 strong component for v1 and v2, followed by explicit equations to compute v3 and der(x). Then, there is another 2x2 strong component for v4 and v5, followed by an explicit equation for v6, which however are not necessary to compute the derivatives.

The code to computeF(z, der(z), t) is thus generated as follows:

F1 := v1 + v2 - u;
F2  := v1 - v2 + u;
v3 := v1 + v2;
F3 := der(x) + x - v3;

in this case, w = [v1 v2]'. Note that IDA will solve a 3x3 system in v1, v2, and x, and doesn't see v3, v4, and v5, nor does not take care of computing them explicitly. The 3x3 symbolic Jacobian for this system will be generated with reference to the first three variables x, v1, v2 and to the three residuals F1, F2, F3 only.

IDA needs to compute F(z, der(z), t) multiple times for each step, until it reaches convergence. It is thus good for performance that F has the minimum possible dimension, 3 in this case. Once it has converged to a solution for x, v1, v2, the code to compute the remaining variables:

v3 := v1 + v2;
solve: v4 + v5 - x = 0
       v4 - v5 + x = 0
for v4, v5
v6 := v4 + v5;

will only be executed once.

Note that if this method is used, purely algebraic systems won't actually use IDA, because there is no need to compute the "next state"; the equations will be solved by the embedded Newton solver of the OMC runtime and the two Jacobians will be empty. For systems that can be fully causalized with explicit assignments, IDA will only see the differential equation residuals, not the algebraic equation residuals.

One the one hand, this method is in general more efficient than throwing the raw DAEs to IDA, so if one seeks maximum efficiency, this is the way to go. On the other hand, grid models are characterized by a very large strong component, containing the network equations, so the advantage of going this way is less dramatic than with systems that have many variables but only a few states and a few implicit equations (say, 3D robot models).

I think that in the medium term it would be advisable that Dynawo can exploit this smarter way of solving the DAEs. In the shorter term, it may be useful if we could skip the causalization step and basically just pass the raw DAEs (with symbolically computed Jacobians) to IDA. As RTE's problems are index-1, this should work out of the box. This was actually Willi Braun's first implementation of daeMode, which was then superseded by the final, more efficient one. For a while, we used to be able to use both methods (see the test results in DOI:10.3384/ecp17132557), then the flags were refactored and eventually the raw DAEs method became no longer available, --daeMode currently enables the causalization-based method only.

@kabdelhak, @phannebohm could add some flag, e.g. -d=rawDaeMode to skip the causalization step and just pass the raw residual equations to IDA?

[casella]

One more point to the discussion: the C code that Dynawo needs from OMC at the moment does not correspond to a whole system (that will come when the new array-preserving back-end is ready to handle it), but only to an individual component model with causal (input-output) and a-causal (electrical connectors) interfaces. If this model is a controller model with inputs and outputs, it can be usually causalized completely in explicit form, because control diagrams are not expected to contain algebraic loops and can be computed sequentially. But if the model has an a-causal interface, nothing can be assumed on causality at the boundary, so I guess OMC would see it as a single strong component. In that case, you would get full Jacobians out of the box without further intervention.

I'm not really sure about that, I should understand better what boundary conditions are applied to the a-causal component ports when you generate the C code for Dynawo. I guess you are probably connecting them to some dummy component containing a generic f(port.v, port.i) = 0 equation that you then discard when using the code in Dynawo, but I'm not sure. In that case, no causalization will be possible, so you should get the full Jacobian already. Except, of course, for equations that will end up in the "leaves" of the dependency diagrams, such as, e.g., the computation of phase angles at the ports. But I guess you don't really care about them for your grid simulations, so it could already be OK.

@marcochiaramello, @rosiereflo, @joyelfeghali can you comment on that?