Simulating a pure algebraic model in daeMode leads to a wrong solution

[Francesco Casella]

The DistributionSystemModelica_N_X_M_X models of the ScalableTestSuite run with daeMode = new badly fail the verification phase, see e.g. ​this report. Analysis of the results reveal that despite V_source.v correctly ramping up to 600 V, all voltages and currents in the distribution networks turn out to be zero.

You can replicate the issue with a much simpler model:

model Test
  extends ScalableTestSuite.Electrical.DistributionSystemDC.ScaledExperiments.DistributionSystemModelica_N_10_M_10
   (N = 1, M = 1);

annotation(
  experiment(StopTime = 1, Interval = 1e-3),
  uses(ScalableTestSuite(version = "1.11.4")),
  __OpenModelica_commandLineOptions = "-d=optdaedump --daeMode=new");
end Test;

which feature the series connection of the voltage source V_source, the resistor primary, the resistor secondary, and the resistor load, which is eventually connected to ground. As soon as the voltage source is ramped up, all voltages are expected to grow and current to flow in the series connection. Instead, in the simulation results, primary[1].v is identically zero, so no current flows, and then everything else is zero accordingly.

The last step of optdaedump outputs the following equations

Variables (9) 
======================================== 
1: V_source.v:VARIABLE(unit = "V" )  "Voltage drop between the two pins (= p.v - n.v)".Test, .Modelica.Electrical.Analog.Sources.RampVoltage$V_source type: Real  
2: $DAEres0:DAE_RESIDUAL_VAR()  type: Real  
3: $DAEres1:DAE_RESIDUAL_VAR()  type: Real  
4: $DAEres2:DAE_RESIDUAL_VAR()  type: Real  
5: $DAEres3:DAE_RESIDUAL_VAR()  type: Real  
6: $DAEres4:DAE_RESIDUAL_VAR()  type: Real  
7: primary[1].LossPower:VARIABLE(unit = "W" )  "Loss power leaving component via heatPort".Test, .Modelica.Electrical.Analog.Basic.Resistor$primary type: Real  [1] 
8: load[1,1].LossPower:VARIABLE(unit = "W" )  "Loss power leaving component via heatPort".Test, .Modelica.Electrical.Analog.Basic.Resistor$load type: Real  [1,1] 
9: secondary[1,1].LossPower:VARIABLE(unit = "W" )  "Loss power leaving component via heatPort".Test, .Modelica.Electrical.Analog.Basic.Resistor$secondary type: Real  [1,1] 
 
 
Equations (9, 9) 
======================================== 
1/1 (1): V_source.v := V_source.signalSource.offset + (if time < V_source.signalSource.startTime then 0.0 else if time < V_source.signalSource.startTime + V_source.signalSource.duration then (time - V_source.signalSource.startTime) * V_source.signalSource.height / V_source.signalSource.duration else V_source.signalSource.height)   [auxiliary |1|1|0|1|] 
2/2 (1): $DAEres0 := load[1,1].R_actual * secondary[1,1].i - load[1,1].v   [dynamic |1|0|0|1|] 
3/3 (1): $DAEres1 := secondary[1,1].R_actual * secondary[1,1].i - secondary[1,1].v   [dynamic |1|0|0|1|] 
4/4 (1): $DAEres2 := primary[1].n.v + (-load[1,1].v) - secondary[1,1].v   [dynamic |1|0|0|1|] 
5/5 (1): $DAEres3 := V_source.v + (-primary[1].n.v) - primary[1].v   [dynamic |1|0|0|1|] 
6/6 (1): $DAEres4 := primary[1].R_actual * secondary[1,1].i - primary[1].v   [dynamic |1|0|0|1|] 
7/7 (1): primary[1].LossPower := primary[1].v * secondary[1,1].i   [auxiliary |0|1|0|1|] 
8/8 (1): load[1,1].LossPower := load[1,1].v * secondary[1,1].i   [auxiliary |0|1|0|1|] 
9/9 (1): secondary[1,1].LossPower := secondary[1,1].v * secondary[1,1].i   [auxiliary |0|1|0|1|] 

Equation 5/5 is correct, and states that the voltage across the source V_source.v should be equal to the sum of the voltage across the primary resistor primary[1].v plus the voltage on its negative pin primary[1].n.v. If you check the simulation results, though, the first voltage ramps up to 600, while the other two remain zero.

I guess something goes wrong in the code generation. Test_16dae.c contains the following lines

/*
 equation index: 58
 type: SIMPLE_ASSIGN
 $DAEres3 = V_source.v + (-primary[1].n.v) - primary[1].v
 */
void Test_eqFunction_58(DATA *data, threadData_t *threadData)
{
  TRACE_PUSH
  const int equationIndexes[2] = {1,58};
  $P$DAEres3 = data->localData[0]->realVars[0] /* V_source.v variable */ + (-data->localData[0]->realVars[7] /* primary[1].n.v variable */) - data->localData[0]->realVars[8] /* primary[1].v variable */;
  TRACE_POP
}

which seems correct looking at the comments.

However, the Test_info.json file contains this chunk:

"variables":{
"V_source.v":{"comment":"Voltage drop between the two pins (= p.v - n.v)","kind":"variable","type":"Real","unit":"V","displayUnit":"","source":{"info":{"file":"C:/OpenModelica1.13.0-dev-64bit/lib/omlibrary/Modelica 3.2.2/Electrical/Analog/Interfaces.mo","lineStart":173,"lineEnd":173,"colStart":5,"colEnd":67},"within":["Real"],"instance":"V_source","typeLst":["Test","Modelica.Electrical.Analog.Sources.RampVoltage$V_source"]}},
"load[1,1].LossPower":{"comment":"Loss power leaving component via heatPort","kind":"variable","type":"Real","unit":"W","displayUnit":"","source":{"info":{"file":"C:/OpenModelica1.13.0-dev-64bit/lib/omlibrary/Modelica 3.2.2/Electrical/Analog/Interfaces.mo","lineStart":312,"lineEnd":312,"colStart":5,"colEnd":67},"within":["Real"],"instance":"load[1,1]","typeLst":["Test","Modelica.Electrical.Analog.Basic.Resistor$load"]}},
"load[1,1].R_actual":{"comment":"Actual resistance = R*(1 + alpha*(T_heatPort - T_ref))","kind":"variable","type":"Real","unit":"Ohm","displayUnit":"","source":{"info":{"file":"C:/OpenModelica1.13.0-dev-64bit/lib/omlibrary/Modelica 3.2.2/Electrical/Analog/Basic.mo","lineStart":62,"lineEnd":63,"colStart":5,"colEnd":63},"within":["Real"],"instance":"load[1,1]","typeLst":["Test","Modelica.Electrical.Analog.Basic.Resistor$load"],"operations":[{"op":"info","display":"solved","data":["load[1,1].R_actual = load[1,1].R * (1.0 + load[1,1].alpha * (load[1,1].T - load[1,1].T_ref))"]},{"op":"chain","display":"substitution","data":["load[1,1].R * (1.0 + load[1,1].alpha * (load[1,1].T_heatPort - load[1,1].T_ref))","load[1,1].R * (1.0 + load[1,1].alpha * (load[1,1].T - load[1,1].T_ref))"]},{"op":"before-after","display":"flattening","data":["R_actual = R * (1 + alpha * (T_heatPort - T_ref));","load[1,1].R_actual = load[1,1].R * (1.0 + load[1,1].alpha * (load[1,1].T_heatPort - load[1,1].T_ref));"]}]}},
"load[1,1].v":{"comment":"Voltage drop between the two pins (= p.v - n.v)","kind":"variable","type":"Real","unit":"V","displayUnit":"","source":{"info":{"file":"C:/OpenModelica1.13.0-dev-64bit/lib/omlibrary/Modelica 3.2.2/Electrical/Analog/Interfaces.mo","lineStart":173,"lineEnd":173,"colStart":5,"colEnd":67},"within":["Real"],"instance":"load[1,1]","typeLst":["Test","Modelica.Electrical.Analog.Basic.Resistor$load"]}},
"primary[1].LossPower":{"comment":"Loss power leaving component via heatPort","kind":"variable","type":"Real","unit":"W","displayUnit":"","source":{"info":{"file":"C:/OpenModelica1.13.0-dev-64bit/lib/omlibrary/Modelica 3.2.2/Electrical/Analog/Interfaces.mo","lineStart":312,"lineEnd":312,"colStart":5,"colEnd":67},"within":["Real"],"instance":"primary[1]","typeLst":["Test","Modelica.Electrical.Analog.Basic.Resistor$primary"]}},
"primary[1].R_actual":{"comment":"Actual resistance = R*(1 + alpha*(T_heatPort - T_ref))","kind":"variable","type":"Real","unit":"Ohm","displayUnit":"","source":{"info":{"file":"C:/OpenModelica1.13.0-dev-64bit/lib/omlibrary/Modelica 3.2.2/Electrical/Analog/Basic.mo","lineStart":62,"lineEnd":63,"colStart":5,"colEnd":63},"within":["Real"],"instance":"primary[1]","typeLst":["Test","Modelica.Electrical.Analog.Basic.Resistor$primary"],"operations":[{"op":"info","display":"solved","data":["primary[1].R_actual = primary[1].R * (1.0 + primary[1].alpha * (primary[1].T - primary[1].T_ref))"]},{"op":"chain","display":"substitution","data":["primary[1].R * (1.0 + primary[1].alpha * (primary[1].T_heatPort - primary[1].T_ref))","primary[1].R * (1.0 + primary[1].alpha * (primary[1].T - primary[1].T_ref))"]},{"op":"before-after","display":"flattening","data":["R_actual = R * (1 + alpha * (T_heatPort - T_ref));","primary[1].R_actual = primary[1].R * (1.0 + primary[1].alpha * (primary[1].T_heatPort - primary[1].T_ref));"]}]}},
"primary[1].n.v":{"comment":"Potential at the pin","kind":"variable","type":"Real","unit":"V","displayUnit":"","source":{"info":{"file":"C:/OpenModelica1.13.0-dev-64bit/lib/omlibrary/Modelica 3.2.2/Electrical/Analog/Interfaces.mo","lineStart":87,"lineEnd":92,"colStart":5,"colEnd":61},"within":["Real"],"instance":"primary[1]n","typeLst":["Test","Modelica.Electrical.Analog.Basic.Resistor$primary","Modelica.Electrical.Analog.Interfaces.NegativePin"]}},
"primary[1].v":{"comment":"Voltage drop between the two pins (= p.v - n.v)","kind":"variable","type":"Real","unit":"V","displayUnit":"","source":{"info":{"file":"C:/OpenModelica1.13.0-dev-64bit/lib/omlibrary/Modelica 3.2.2/Electrical/Analog/Interfaces.mo","lineStart":173,"lineEnd":173,"colStart":5,"colEnd":67},"within":["Real"],"instance":"primary[1]","typeLst":["Test","Modelica.Electrical.Analog.Basic.Resistor$primary"]}},
"secondary[1,1].LossPower":{"comment":"Loss power leaving component via heatPort","kind":"variable","type":"Real","unit":"W","displayUnit":"","source":{"info":{"file":"C:/OpenModelica1.13.0-dev-64bit/lib/omlibrary/Modelica 3.2.2/Electrical/Analog/Interfaces.mo","lineStart":312,"lineEnd":312,"colStart":5,"colEnd":67},"within":["Real"],"instance":"secondary[1,1]","typeLst":["Test","Modelica.Electrical.Analog.Basic.Resistor$secondary"]}},

and if you count the variables from zero (as usual with C records) it seems to me that
primary[1].n.v and primary[1].v are number 6 and 7, not 7 and 8. I guess that you are mapping Modelica arrays into C arrays without subtracting 1 somewhere in the daeMode=new code generation functions, which also explains why this problem only shows up with vectorized ScalableTestSuite test cases.

Please fix this ASAP, before the 1.13.0 release.

Also, I have to admit I do not fully understand the optdaedump output. The DAE residuals are listed among the unknowns, but they are in fact not unknown at all, as they should be zero at the solution. Conversely, the real unknowns of the corresponding strong component of the BLT, namely primary[1].v, primary[1].n.v, secondary[1,1].v, and secondary[1,1].i do not show up anywhere, which is quite puzzling. Unless I am missing something, please make sure the correct unknowns are listed.

[Willi Braun]

The issue here is not related to any indexing issue, it seems that we just disable the integrator on this model, for some buggy reasons. Should be fixed easily.

[Francesco Casella]

Replying to wbraun:

The issue here is not related to any indexing issue, it seems that we just disable the integrator on this model, for some buggy reasons. Should be fixed easily.

Good!

For my information, I'd nevertheless like to understand why the indeces in the C-code do not match with the json file. Is some other information source used to determine the matching between variable names and data indeces?

[Willi Braun]

Replying to casella:

Replying to wbraun:

The issue here is not related to any indexing issue, it seems that we just disable the integrator on this model, for some buggy reasons. Should be fixed easily.

Good!

Fixed in ​PR2740.

For my information, I'd nevertheless like to understand why the indeces in the C-code do not match with the json file. Is some other information source used to determine the matching between variable names and data indeces?

Yes, you are right the indeces in the json file are wrong, as far as I see we miss there for some reasons a variable: ground[1,1].p.v. This needs some more investigation and probably an other ticket.

[Willi Braun]

Replying to wbraun:

Replying to casella:

Replying to wbraun:

The issue here is not related to any indexing issue, it seems that we just disable the integrator on this model, for some buggy reasons. Should be fixed easily.

Good!

Fixed in ​PR2740.

For my information, I'd nevertheless like to understand why the indeces in the C-code do not match with the json file. Is some other information source used to determine the matching between variable names and data indeces?

Yes, you are right the indeces in the json file are wrong, as far as I see we miss there for some reasons a variable: ground[1,1].p.v. This needs some more investigation and probably an other ticket.

The ticket:3789 reports already the issue with the json file.

[Francesco Casella]

Replying to wbraun:

Thanks for the quick fix!

The ticket:3789 reports already the issue with the json file.

I see. I guess this is not critical. Is there any other place I can look at to see who's who in the data structure?

[Willi Braun]

Replying to casella:

Replying to wbraun:

Thanks for the quick fix!

The ticket:3789 reports already the issue with the json file.

I see. I guess this is not critical. Is there any other place I can look at to see who's who in the data structure?

Basically I'm using mostly the dumps dumpSimCode and dumpindxdae.

[Francesco Casella]

Replying to wbraun:

Basically I'm using mostly the dumps dumpSimCode and dumpindxdae.

Thanks for the hint!

I just started the Hudson job to re-run the daeMode tests, let's see if everything works out properly.