Backend fails with discrete-time first order model

[Rüdiger Franke]

OpenModelica works with lots of discrete-time models, including also more sophisticated ones than here, see e.g. 4 states with implicit Euler:
​https://github.com/OpenModelica/OpenModelica-testsuite/blob/master/openmodelica/cppruntime/fmu/modelExchange/2.0/testCSTR.mos

However, a simple model with one linear state fails. The below model samples a FirstOrder model. -d=backenddaeinfo says:

[7] 18:26:42 Symbolisch Meldung
Torn system details for strict tearing set:
 * Linear torn systems: 1 {(1,100.0%) 1}
 * Non-linear torn systems: 0

-d=dumpSimCode shows the following equation system -- without Jacobian and with index: -1 !?

10:  (LINEAR) index:0 jacobian: false 
        variables: 
index: -1: firstOrder.y (no alias)  initial: firstOrder.y_start no arrCref index:() []  
        b-vector: 
        8: $DER.firstOrder.y=(firstOrder.k * firstOrder.u - firstOrder.y) / firstOrder.T [Real] 
        9: interval() * (if firstTick() then 0.0 else $DER.firstOrder.y) + previous(firstOrder.y) - firstOrder.y (RESIDUAL) 

The C runtime reports "Failed to solve linear system of equations (no. 10) at time 0.000000, system is singular for U[2, 2].".
The Cpp runtime removes the equation during code generation.

model DiscreteTimeLoop
  Modelica.Blocks.Continuous.FirstOrder firstOrder(T = 1)  annotation(
    Placement(visible = true, transformation(origin = {18, 10}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
  Modelica_Synchronous.RealSignals.Sampler.SampleClocked sample1 annotation(
    Placement(visible = true, transformation(origin = {-24, 10}, extent = {{-6, -6}, {6, 6}}, rotation = 0)));
  Modelica_Synchronous.ClockSignals.Clocks.PeriodicRealClock periodicClock1(period = 0.1)  annotation(
    Placement(visible = true, transformation(origin = {-54, -44}, extent = {{-6, -6}, {6, 6}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealInput u(start = 1)  annotation(
    Placement(visible = true, transformation(origin = {-82, 10}, extent = {{-8, -8}, {8, 8}}, rotation = 0), iconTransformation(origin = {-72, 10}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
equation
  connect(u, sample1.u) annotation(
    Line(points = {{-82, 10}, {-32, 10}, {-32, 10}, {-32, 10}}, color = {0, 0, 127}));
  connect(sample1.y, firstOrder.u) annotation(
    Line(points = {{-18, 10}, {6, 10}, {6, 10}, {6, 10}}, color = {0, 0, 127}));
  connect(periodicClock1.y, sample1.clock) annotation(
    Line(points = {{-48, -44}, {-24, -44}, {-24, 2}, {-24, 2}}, color = {175, 175, 175}));

annotation(
    uses(Modelica(version = "3.2.3"), Modelica_Synchronous(version = "0.93.0")));
end DiscreteTimeLoop;

[Francesco Casella]

@rfranke, 1.16.1 has been released already. If this is a critical fix for you, please assign it to 1.16.2

[Rüdiger Franke]

It requires workarounds in lot of our models, not really helping to motivate the use of OpenModelica.
We are speaking of a linear equation with one unknown variable here!?

[Martin Sjölund]

A linear equation with one unknown variable is simply called an equation. These are 2x2 systems.
The model works fine on the C runtime using the 1.16.0 docker image though, so I'm not sure what the actual error is.
Is this a problem on Windows or Linux? Does it occur only on a particular CPU?

[Rüdiger Franke]

Did you also check the result? Constant zero for firstOrder.y or a staircase?
Here is the equation from the opener again. There is no 2x2 system, just one unknown variable firstOrder.y. $DER.firstOrder.y can easily be eliminated.

10:  (LINEAR) index:0 jacobian: false 
        variables: 
index: -1: firstOrder.y (no alias)  initial: firstOrder.y_start no arrCref index:() []  
        b-vector: 
        8: $DER.firstOrder.y=(firstOrder.k * firstOrder.u - firstOrder.y) / firstOrder.T [Real] 
        9: interval() * (if firstTick() then 0.0 else $DER.firstOrder.y) + previous(firstOrder.y) - firstOrder.y (RESIDUAL) 

[Francesco Casella]

Testing with v1.17.0-dev-220-g67bec2f83f, I get the same problem with the C-runtime

As you can see from the debugger, the system has two equations, the unknowns are

firstOrder.y
der(firstOrder.y)

The system is currently solved by tearing in the backend:

(torn) der(firstOrder.y) := (firstOrder.k*firstOrder.u - firstOrder.y)/firstOrder.T
(residual) interval()*(if firstTick() then 0.0 else $DER.firstOrder.y) +
           previous(firstOrder.y) - firstOrder.y = 0

When firstTick() is true the residual is zero for

firstOrder.y = previous(firstOrder.y)

then the torn equation computes

der(firstOrder.y) := (firstOrder.k*firstOrder.u - firstOrder.y)/firstOrder.T

This works fine, if you set StopTime = 0 the simulation is successful. At the next tick the variable $DER.firstOrder.y comes into play, which is quite confusing, but that's already the topic of #6150. I understand that variable actually means der(firstOrder.y), so the residual equation would eventually be equivalent to

0.1*firstOrder.u - 0.1*firstOrder.y + previous(firstOrder.y) - firstOrder.y = 0

which is not singular. The solution should be

firstOrder.y = (0.1*firstOrder + previous(firstOrder.y))/ 1.1

Maybe this is related to #6196?

@AnHeuermann, could you please have a look and check what is actually written in the A and b matrices for the linear solver?

Please back-port any fix that could solve this issue to the maintenance branch as well.

BTW, I'm still unsure what we'll do with linear tearing (see #6196). In any case, if the system is torn to just one linear equation, there is no need to solve that using dgesv(), we could solve that directly as

x := b/A

[Vitalij Ruge]

Replying to rfranke:

Did you also check the result? Constant zero for firstOrder.y or a staircase?
Here is the equation from the opener again. There is no 2x2 system,

For the sake of completeness, the workaround --maxSizeSolveLinearSystem=2 work for this small model with c/cpp. Of course the risk is high, the workaround fail for real problem.

[Karim Adbdelhak]

The following .mos file works for me on the latest master (OMCompiler v1.17.0-dev.241+gf877feaade) for both C and C++:

loadFile("DiscreteTimeLoop.mo"); getErrorString();

setCommandLineOptions("--std=3.3 -d=newInst");
simulate(DiscreteTimeLoop); getErrorString();

This is the reported algebraic loop using -d=dumpLoops which seems perfectly fine. I guess this has been fixed by a fairly recent change i did in ​PR6985 or ​PR6988

########################################
                                                                                                     
 Algbraic Loops (Simulation):         
                                                                                                     
########################################                                                             
                                                                                                     
system 1                                                                                             
========================================
                                                                                                     
################################################################################
 dumpLoops: SORTED COMPONENT 
################################################################################
                                                  
torn linear Equationsystem:
                                                  
internal vars (1)                                                                                    
1: $DER.firstOrder.y:VARIABLE()  type: Real
                                                                                                     
residual vars (1)                                                                                    
1: firstOrder.y:CLOCKED_STATE(flow=false start = firstOrder.y_start fixed = true )  "Connector of Real output signal" type: Real
                                                                                                     
internal equations (1)                                                                               
1/1 (1): $DER.firstOrder.y = (firstOrder.k * firstOrder.u - firstOrder.y) / firstOrder.T   [dynamic |0|0|0|0|]
                                                                                                     
residual equations (1)
1/1 (1): (firstOrder.y - previous(firstOrder.y)) / interval() = if firstTick() then 0.0 else $DER.firstOrder.y   [dynamic |0|0|0|0|]

I will close this ticket, please reopen if there are still issues!

[Karim Adbdelhak]

Note that it also works with --tearingMethod=minimalTearing, which should reflect the original problem:

########################################
                                                                                                     
 Algbraic Loops (Simulation):         
                                                                                                     
########################################                                                             
                                                                                                     
system 1                                                                                             
========================================
                                                                                                     
################################################################################
 dumpLoops: SORTED COMPONENT 
################################################################################
                                                  
torn linear Equationsystem:
                                                  
internal vars (0)                                                                                    
                                                  
residual vars (2)                                                                                    
1: $DER.firstOrder.y:VARIABLE()  type: Real                                                          
2: firstOrder.y:CLOCKED_STATE(flow=false start = firstOrder.y_start fixed = true )  "Connector of Real output signal" type: Real
                                                                                                     
internal equations (0)                                                                               
                                                                                                     
residual equations (2)                                                                               
1/1 (1): (firstOrder.y - previous(firstOrder.y)) / interval() = if firstTick() then 0.0 else $DER.firstOrder.y   [dynamic |0|0|0|0|]
2/2 (1): $DER.firstOrder.y = (firstOrder.k * firstOrder.u - firstOrder.y) / firstOrder.T   [dynamic |0|0|0|0|]

[Francesco Casella]

@Karim, @rfranke needs this fix in a stable release, see comment:4, so it should be back-ported to the 1.16 maintenance branch, and then released as 1.16.2. Can you please check that patching your two PR's onto the maintenance branch resolve the issue and then push them there as well?

Thanks!

[Karim Adbdelhak]

I tested it with the v1.16 maintenance branch (OMCompiler v1.16.1-v1.16.1.5+gc01d479c53) and it reports the same. Maybe it got fixed before? I cannot reproduce the problem.

[Rüdiger Franke]

Strange: the C runtime generates no error when running your omc script. This confirms the comment by @sjoelund.se above.

The C runtime raises an error when running from OMEdit though!?

Failed to solve linear system of equations (no. 10) at time 0.000000, system is singular for U[2, 2].
The default linear solver fails, the fallback solver with total pivoting is started at time 0.000000. That might raise performance issues, for more information use -lv LOG_LS.
Matrix singular!

But also when running from the command line, the result of C and C++ runtime is wrong (firstOrder.y constant 1 instead of gradually raising from 0 to 0.67).

[Francesco Casella]

So: runtime fails with OMEdit giving singular system, and gives wrong result from a script. Interesting.

@Karim, any clue what's going on?

[Francesco Casella]

Ping!

Sorry to push you, can we fix this in 1.16.2?

[Andreas Heuermann]

But also when running from the command line, the result of C and C++ runtime is wrong (firstOrder.y constant 1 instead of gradually raising from 0 to 0.67).

So I tested with the following script:

setCommandLineOptions("-d=newInst --simCodeTarget=Cpp --std=3.3"); getErrorString();
loadFile("DiscreteTimeLoop.mo"); getErrorString();
simulate(DiscreteTimeLoop); getErrorString();

This will simulate, but the result is wrong. firstOrder.y is constant 0.

---

So I tested with --simCodeTarget=C and got:

true
""
true
"Notification: Automatically loaded package Modelica 3.2.3 due to uses annotation.
Notification: Automatically loaded package Complex 3.2.3 due to uses annotation.
Notification: Automatically loaded package ModelicaServices 3.2.3 due to uses annotation.
Notification: Automatically loaded package Modelica_Synchronous 0.93.0 due to uses annotation.
"
record SimulationResult
    resultFile = "/home/andreas/workspace/Testitesttest/tickets/ticket6254/DiscreteTimeLoop_res.mat",
    simulationOptions = "startTime = 0.0, stopTime = 1.0, numberOfIntervals = 500, tolerance = 1e-06, method = 'dassl', fileNamePrefix = 'DiscreteTimeLoop', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags = ''",
    messages = "LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
stdout            | warning | The default linear solver fails, the fallback solver with total pivoting is started at time 0.000000. That might raise performance issues, for more information use -lv LOG_LS.
LOG_SUCCESS       | info    | The simulation finished successfully.
",
    timeFrontend = 0.6425105,
    timeBackend = 0.0525949,
    timeSimCode = 0.0118169,
    timeTemplates = 0.0107241,
    timeCompile = 2.4457127,
    timeSimulation = 0.0308921,
    timeTotal = 3.1944813
end SimulationResult;
""
"/home/andreas/workspace/Testitesttest/tickets/ticket6254/DiscreteTimeLoop.fmu"
""

Here the result for firstOrder.y is constant 0 as well.
Since both Cpp and C doesn't seem to work I'll concentrate on the C runtime for now.

---

So looking at the generated C code we have a strange "empty" loop:

/* linear systems */
OMC_DISABLE_OPT
void setLinearMatrixA10(void *inData, threadData_t *threadData, void *systemData)
{
  const int equationIndexes[2] = {1,10};
  DATA* data = (DATA*) inData;
  LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData;
}
OMC_DISABLE_OPT
void setLinearVectorb10(void *inData, threadData_t *threadData, void *systemData)
{
  const int equationIndexes[2] = {1,10};
  DATA* data = (DATA*) inData;
  LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData;
}
OMC_DISABLE_OPT
void initializeStaticLSData10(void *inData, threadData_t *threadData, void *systemData)
{
  DATA* data = (DATA*) inData;
  LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData;
  int i=0;
  /* static ls data for firstOrder.y */
  linearSystemData->nominal[i] = data->modelData->realVarsData[3].attribute /* firstOrder.y */.nominal;
  linearSystemData->min[i]     = data->modelData->realVarsData[3].attribute /* firstOrder.y */.min;
  linearSystemData->max[i++]   = data->modelData->realVarsData[3].attribute /* firstOrder.y */.max;
}

but -d=dumpLoops tells us we should have a loop generated from:

################################################################################
 dumpLoops: SORTED COMPONENT 
################################################################################

torn linear Equationsystem:

internal vars (1)
1: $DER.firstOrder.y:VARIABLE()  type: Real

residual vars (1)
1: firstOrder.y:CLOCKED_STATE(flow=false start = firstOrder.y_start fixed = true )  "Connector of Real output signal" type: Real

internal equations (1)
1/1 (1): $DER.firstOrder.y = (firstOrder.k * firstOrder.u - firstOrder.y) / firstOrder.T   [dynamic |0|0|0|0|]

residual equations (1)
1/1 (1): (firstOrder.y - previous(firstOrder.y)) / interval() = if firstTick() then 0.0 else $DER.firstOrder.y   [dynamic |0|0|0|0|]

And as @rfranke stated in the ticket description the simCodeDump is looking fishy:

*********************
* SimCode Equations *
*********************


allEquations:
========================================

========================================



odeEquations (0 systems):
========================================
========================================



algebraicEquations (0 systems):
========================================
========================================


clockPartitions (1 systems):

BaseClock:
========================================
Clock(RCONST 0.1
)========================================

SubPartition Vars:
========================================
index:-1: $DER.firstOrder.y (no alias)  initial: 	no arrCref index:() [] 
index:-1: firstOrder.u (no alias)  initial: 	no arrCref index:() [] 
index:-1: firstOrder.y (no alias)  initial: firstOrder.y_start	no arrCref index:() [] (PREVIOUS)
index:-1: sample1.y (no alias)  initial: 	no arrCref index:() [] 
partition equations:
========================================
11: $CLKPRE.firstOrder.y=firstOrder.y [Real]
6: sample1.y=$getPart(u) [Real]
7: firstOrder.u=sample1.y [Real]
10:  (LINEAR) index:0 jacobian: false
	variables:
index:-1: firstOrder.y (no alias)  initial: firstOrder.y_start	no arrCref index:() []
	b-vector:
	8: $DER.firstOrder.y=(firstOrder.k * firstOrder.u - firstOrder.y) / firstOrder.T [Real]
	9: interval() * (if firstTick() then 0.0 else $DER.firstOrder.y) + previous(firstOrder.y) - firstOrder.y (RESIDUAL)

So my guess: We are erroring during code generation but don't report the error.

[Andreas Heuermann]

So when I dump the hash table crefToSimVarHT, that is created in createSimCode and which is passed to the templates, we get:

HashTable:
{$CLKPRE.firstOrder.y,{#SimVar(index=0,name=$CLKPRE.firstOrder.y)#}}
{$DER.firstOrder.y,{#SimVar(index=1,name=$DER.firstOrder.y)#}}
{firstOrder.u,{#SimVar(index=2,name=firstOrder.u)#}}
{firstOrder.y,{#SimVar(index=3,name=firstOrder.y)#}}
{sample1.y,{#SimVar(index=4,name=sample1.y)#}}
{u,{#SimVar(index=5,name=u)#}}
{firstOrder.T,{#SimVar(index=0,name=firstOrder.T)#}}
{firstOrder.k,{#SimVar(index=1,name=firstOrder.k)#}}
{firstOrder.y_start,{#SimVar(index=2,name=firstOrder.y_start)#}}
{periodicClock1.period,{#SimVar(index=3,name=periodicClock1.period)#}}
{firstOrder.initType,{#SimVar(index=0,name=firstOrder.initType)#}}
{periodicClock1.useSolver,{#SimVar(index=0,name=periodicClock1.useSolver)#}}
{sample1.u,{#SimVar(index=0,name=sample1.u)#}}
{periodicClock1.solverMethod,{#SimVar(index=0,name=periodicClock1.solverMethod)#}}

This dump is done directly before the output from -d=dumpSimCode are displayed. So the hash table itself is looking good, but it seems we still can't find some variables (all with index=-1).
If I dump the crefToClockIndexHT as well we get:

HashTable:
{sample1.y,{1}}
{$CLKPRE.sample1.y,{1}}
{firstOrder.y,{1}}
{$CLKPRE.firstOrder.y,{1}}
{firstOrder.u,{1}}
{$CLKPRE.firstOrder.u,{1}}
{$DER.firstOrder.y,{1}}
{$CLKPRE.$DER.firstOrder.y,{1}}

So the next step on debugging this would be to check why in function dumpSimCodeDebug we get index=-1.

[Francesco Casella]

@AnHeuermann will work on that this week, if he can find a solution, we'll put it in 1.16.2, otherwise we'll try to fix that for 1.17.0.

[Andreas Heuermann]

@kabdelhak and myself managed to find the cause for the problem.
The function translateClockedEquations will translate equations from the Backend structure to the SimCode structure. On of it's sub-functions createEquationsForSystem1 is missing a case for BackendDAE.TORNSYSTEM() to handle this torn equation system.

We think we can add the missing case by using createTornSystem. It will probably need some adaptions to run, maybe some code generation is missing as well. But now we have a point to start working on it.

I'm working on it on my branch: ​https://github.com/AnHeuermann/OpenModelica/tree/ticket6254-clockedTornSystems

[Francesco Casella]

Great, thanks!

[Andreas Heuermann]

Today I'm not so optimistic that we actually found the root cause.
It seems the reported index=-1 is normal for torn systems, the simCodeDump just isn't very good.

Comparing to some random example with a torn linear system I noticed, that the simCodeDump is reporting index=-1 for the variables in the linear system as well.
One big difference though is that the working example has a Jacobian, while the clocked system has not.

Maybe the lack of a Jacobian in SimCode is the reason why we don't generate code for the system.

[Andreas Heuermann]

It seems that the lack of the Jacobina is precisely the problem here. In the code generation we loop over list<tuple<Integer, Integer, SimEqSystem>> simJac; to generate the A matrix of a linear system A*x=b, but for our example this is an empty list.

---

In the meantime I also tried to see what will happen for a non-torn linear system. With that we will generate something, but it will faile during compilation of the C code:

DiscreteTimeLoop_03lsy.c:17:89: error: use of undeclared identifier 'clockIndex'
  linearSystemData->setAElement(1, 0, DIVISION_SIM(1.0,data->simulationInfo->clocksData[clockIndex].interval,\"interval()\",equationIndexes), 2, linearSystemData, threadData);
                                                                                        ^
DiscreteTimeLoop_03lsy.c:18:77: error: use of undeclared identifier 'clockIndex'
  linearSystemData->setAElement(1, 1, (-(((data->simulationInfo->clocksData[clockIndex].cnt == 0)?0.0:1.0))), 3, linearSystemData, threadData);
                                                                            ^
DiscreteTimeLoop_03lsy.c:27:152: error: use of undeclared identifier 'clockIndex'
  linearSystemData->setBElement(1, DIVISION_SIM(data->localData[0]->realVars[0] /* previous(firstOrder.y) variable */,data->simulationInfo->clocksData[clockIndex].interval,\"interval()\",equationIndexes), linearSystemData, threadData);
                                                                                                                                                       ^
3 errors generated.
make: *** [<builtin>: DiscreteTimeLoop_03lsy.o] Error 1
make: *** Waiting for unfinished jobs....

The Cpp runtime has a similar error.

---

What we could do is:

• Fix code generation for non-torn linear systems with clocked state variables.
  • This could be a fix for simple models, but more complicated loops with discrete variables will need tearing.
• Enable torn systems to work without symbolic Jacobians.
  • They are not supported at the moment. We would need to rewrite some parts on how we generate linear systems in code generation.

@kabdelhak Any idea how we could do this?
I have to understand clocked systems and how we are generating SimCode structure from BackendStructure before starting to code anything.

It will not be possible to fix this this week, so I'm postponing it to 1.17 for now.

[Francesco Casella]

Retargeted to 1.18.0 because of 1.17.0 timed release.

[Francesco Casella]

Ticket retargeted after milestone closed