Bad treatment of algebraic loops in synchronous models

[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.

[lochel]

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

[wbraun]

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

[wbraun]

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?

[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.

[rfranke]

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

[wbraun]

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.

[wbraun]

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?

[rfranke]

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?

[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.

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,

[wbraun]

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?

[rfranke]

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.

[wbraun]

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?

[lochel]

Replying to wbraun:

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

MCP-0024

[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.

[wbraun]

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.

[casella]

Milestone moved to 1.13.0 due to 1.12.0 already being released.

[wbraun]

@rfranke Can mark this as fixed?

[rfranke]

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?

[wbraun]

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