Bad treatment of algebraic loops in synchronous models

[Rüdiger Franke]

​PR1391 along with ​commit d32ecb make ControlledMixingUnit work by applying index reduction and appropriately implementing ExplicitEuler, respectively.

Some more issues show up when replacing ExplicitEuler with ImplicitEuler, introducing algebraic loops. See the following example:

model ControlledMixingUnit2
  extends Modelica_Synchronous.Examples.Systems.ControlledMixingUnit(
    periodicClock1(solverMethod="ImplicitEuler"));
end ControlledMixingUnit2;

No Jacobians and no ModelStructure are generated!

Moreover, there is a linear torn system of dimension 1 -- shouldn't this be solved analytically?

Last but not least, the C runtime does not generate code for algebraic loops in clocked partitions.

[Lennart Ochel]

Patrick, can you please have a look at the torn system?

[Willi Braun]

Actually they are generated, but currently not processed by SimCodeUtil.

[Willi Braun]

Replying to rfranke:

​PR1391 along with ​commit d32ecb make ControlledMixingUnit work by applying index reduction and appropriately implementing ExplicitEuler, respectively.

Some more issues show up when replacing ExplicitEuler with ImplicitEuler, introducing algebraic loops. See the following example:

model ControlledMixingUnit2
  extends Modelica_Synchronous.Examples.Systems.ControlledMixingUnit(
    periodicClock1(solverMethod="ImplicitEuler"));
end ControlledMixingUnit2;

No Jacobians and no ModelStructure are generated!

Moreover, there is a linear torn system of dimension 1 -- shouldn't this be solved analytically?

Last but not least, the C runtime does not generate code for algebraic loops in clocked partitions.

The Jacobians are now also generated at SimCode phase see ​commit 07d95e.

Now the generated code fails, because of unknown identifier:
error: use of undeclared identifier 'clockIndex'

How it is set in the clock partition equations?
Where is this value stored currently?

[Rüdiger Franke]

clockIndex is generated during code generation -- the index of the clocked partition / equation system. It is used to generate code for operators like interval(), looking up the interval of a given clock.

[Rüdiger Franke]

Btw. you will also need to merge ​PR1391 to treat the high-index problem in the clocked partition.

[Willi Braun]

Replying to rfranke:

Btw. you will also need to merge ​PR1391 to treat the high-index problem in the clocked partition.

My commit is based on that PR.

[Willi Braun]

Replying to rfranke:

clockIndex is generated during code generation -- the index of the clocked partition / equation system. It is used to generate code for operators like interval(), looking up the interval of a given clock.

Okay, is it possible to evaluate this index while compile time?

[Rüdiger Franke]

So far things like interval are treated during code generation, along with event processing. The Cpp runtime generates a variable

  const int clockIndex = <%index%>;

into each function that contains code for a clocked partition -- not yet in the Jacobian code.

Your new ​commit 07d95e loops through the clocked partitions. Would it be possible to store the iteration counter as one based clockIndex in the Jacobian?

[Rüdiger Franke]

It appears painful to add something to a Jacibian, because it is defined as tuple, i.e. would require changes at many places. Alternatively the Jacobian vars can be added to the already existing SimCode.crefToClockIndexHT -- e.g. used for clockIndex to modelDescription.xml.

This patch extends your ​commit 07d95e.
It works with an addition to the code generation like: ​Generate Jacobians for clocked partitions in Cpp runtime

SimCodeUtil.mo

@@ -187 +187 @@
   SimCode.HashTableCrefToSimVar crefToSimVarHT;
   SimCodeFunction.MakefileParams makefileParams;
   SimCode.ModelInfo modelInfo;
-  HashTable.HashTable crefToClockIndexHT;
+  HashTable.HashTable crefToClockIndexHT = HashTable.emptyHashTable();
   array<Integer> systemIndexMap;
   list<BackendDAE.EqSystem> clockedSysts, contSysts;
   //list<BackendDAE.Equation> paramAsserts, remEqLst;
@@ -384 +384 @@
     SymbolicJacsNLS := listAppend(SymbolicJacsTemp, SymbolicJacsNLS);
     (allEquations, numberofLinearSys, numberofNonLinearSys, numberofMixedSys, numberOfJacobians, SymbolicJacsTemp) := countandIndexAlgebraicLoops(allEquations, numberofLinearSys, numberofNonLinearSys, numberofMixedSys, numberOfJacobians, {});
     SymbolicJacsNLS := listAppend(SymbolicJacsTemp, SymbolicJacsNLS);
-    (clockedPartitions, numberofLinearSys, numberofNonLinearSys, numberofMixedSys, numberOfJacobians, SymbolicJacsTemp) := countandIndexAlgebraicLoopsClockPartitions(clockedPartitions, numberofLinearSys, numberofNonLinearSys, numberofMixedSys, numberOfJacobians, {});
+    (clockedPartitions, numberofLinearSys, numberofNonLinearSys, numberofMixedSys, numberOfJacobians, SymbolicJacsTemp, crefToClockIndexHT) := countandIndexAlgebraicLoopsClockPartitions(clockedPartitions, numberofLinearSys, numberofNonLinearSys, numberofMixedSys, numberOfJacobians, {}, crefToClockIndexHT);
     SymbolicJacsNLS := listAppend(SymbolicJacsTemp, SymbolicJacsNLS);
 
     // collect symbolic jacobians from state selection
@@ -472 +472 @@
     //BaseHashTable.dumpHashTable(varToArrayIndexMapping);
     //print("END MAPPING\n\n");
 
-    (crefToClockIndexHT, _) := List.fold(inBackendDAE.eqs, collectClockedVars, (HashTable.emptyHashTable(), 1));
+    (crefToClockIndexHT, _) := List.fold(inBackendDAE.eqs, collectClockedVars, (crefToClockIndexHT, 1));
 
     simCode := SimCode.SIMCODE(modelInfo,
                               {}, // Set by the traversal below...
@@ -1200 +1200 @@
   input Integer inMixedSysIndex;
   input Integer inJacobianIndex;
   input list<SimCode.JacobianMatrix> inSymJacs;
+  input HashTable.HashTable inClkHT;
   output list<SimCode.ClockedPartition> outColockPartition = {};
   output Integer outLinearSysIndex = inLinearSysIndex;
   output Integer outNonLinSysIndex = inNonLinSysIndex;
   output Integer outMixedSysIndex = inMixedSysIndex;
   output Integer outJacobianIndex = inJacobianIndex;
   output list<SimCode.JacobianMatrix> outSymJacs = inSymJacs;
+  output HashTable.HashTable outClkHT = inClkHT;
 protected
   SimCode.SubPartition subPartitionTmp;
   SimCode.ClockedPartition clockPartitionTmp;
   list<SimCode.SubPartition> subPartitionsTmp;
   list<SimCode.SimEqSystem> equations;
+  SimCode.JacobianMatrix symJac;
+  list<SimCode.JacobianColumn> columns;
+  list<SimCodeVar.SimVar> simVars;
+  Integer firstJacIdx, clockIndex = 0;
 algorithm
   for clockPartion in inColockPartition loop
     subPartitionsTmp := {};
     for subPartion in clockPartion.subPartitions loop
+      clockIndex := clockIndex + 1;
+      firstJacIdx := listLength(outSymJacs) + 1;
       (equations, outLinearSysIndex, outNonLinSysIndex, outMixedSysIndex, outJacobianIndex, outSymJacs) :=
        countandIndexAlgebraicLoops(subPartion.equations, outLinearSysIndex, outNonLinSysIndex, outMixedSysIndex, outJacobianIndex, outSymJacs);
+      // add jac vars with clockIndex to outClkHT
+      for i in firstJacIdx:listLength(outSymJacs) loop
+        symJac := listGet(outSymJacs, i);
+        (columns, simVars, _, _, _, _, _) := symJac;
+        for simVar in simVars loop
+          outClkHT := BaseHashTable.add((simVar.name, clockIndex), outClkHT);
+        end for;
+        for column in columns loop
+          (_, simVars, _) := column;
+          for simVar in simVars loop
+            outClkHT := BaseHashTable.add((simVar.name, clockIndex), outClkHT);
+          end for;
+        end for;
+      end for;
       subPartitionTmp := SimCode.SUBPARTITION(
         subPartion.vars,
         equations,

[Willi Braun]

Replying to rfranke:

It appears painful to add something to a Jacobian, because it is defined as tuple, i.e. would require changes at many places. Alternatively the Jacobian vars can be added to the already existing SimCode.crefToClockIndexHT -- e.g. used for clockIndex to modelDescription.xml.

True, therefor I started on Friday to rewrite them into records and add the clockIndex to the jacobian record, I'll add this later the day.
Since the clockIndex is a partition index, it's better to store it once in the jacobian, are you prefer your solution with the variables?

[Rüdiger Franke]

A Jacobian record with member clockIndex would of course be preferred. The variables need the clockIndex if they appear in modelDescription.xml. Once the clockIndex is available in the Jacobian, a separate function could transfer that info to crefToClockIndexHT on demand.

[Willi Braun]

In ​PR1439 the jacobian got partition index support.

In the c-runtime we the function "firstTick()" is not found in the jacobians.
@lochel Where is this defined?

[Lennart Ochel]

Replying to wbraun:

In the c-runtime we the function "firstTick()" is not found in the jacobians.
@lochel Where is this defined?

MCP-0024

[Rüdiger Franke]

With ​fixing ticket:4237: Support of clockPartitions in jacobians the example ControlledMixingUnit2 gets clock indices 1, 2 and 3. It should only be 1. Moreover, continuous-time models get a clockIndex 3, which is wrong -- it should be 0 if appearing at all in non-synchronous model code.

[Willi Braun]

Replying to rfranke:

With ​fixing ticket:4237: Support of clockPartitions in jacobians the example ControlledMixingUnit2 gets clock indices 1, 2 and 3. It should only be 1. Moreover, continuous-time models get a clockIndex 3, which is wrong -- it should be 0 if appearing at all in non-synchronous model code.

Okay, should be fixed with ​PR1442.

[Francesco Casella]

Milestone moved to 1.13.0 due to 1.12.0 already being released.

[Willi Braun]

@rfranke Can mark this as fixed?

[Rüdiger Franke]

Yes, this ticket can be closed.

Just opened the new ticket #4783 about a new issue.

Moreover, the FMI ModelStructure is still missing for clocked partitions. Will open a new ticket as well -- don't we already have another ticket for this?

[Willi Braun]

Ah, yes Regarding the discrete structure I've opened #4424 again.