Unknowns with explicitly set start attributes should be preferred as tearing variables

[Francesco Casella]

Consider the model ​Modelica.Electrical.Analog.Examples.OpAmps.SignalGenerator.

The component opAmp1 has a positive feedback connection, with two possible initial equilibrium points: one with v_out = -15, the other one with v_out = 0. This is the result of the solution of a strongly nonlinear system of equations. The following declaration in the model

  Modelica.Electrical.Analog.Basic.Resistor r2(R=R2, i(start=Vps/R2))

leads to the convergence to the desired initial condition (i.e., v_out = -15), which starts a cycle of oscillations, as in the reference solution, as long as r2.i is selected as iteration/tearing variable. Unfortunately, OMC selects a different unknown as tearing variable, which has the default start value of zero; as a consequence, the wrong solution is selected, and no oscillations ensue.

One solution to this issue is to add appropriate start values to all potential iteration variables in the MSL source code, in order to guarantee that the correct solution is found irrespective of the selected iteration variables. However, this is a bit cumbersome.

It would be nicer if the tearing algorithm could give higher priority to r2.i in the selection as tearing variable candidate, since its start value is explicitly modified in the model, while the other potential tearing variables only have the default set by their type.

This feature would make OMC more robust in handling models that were developed with other tools (as in this case), whereby some start values are set in order to steer the convergence of the solver towards the desired solution.

Karim, Andreas, Lennart, do you think this could be implemented without too much effort?

[Lennart Ochel]

This is already implemented in the tearing algorithm, but disabled by default. It can be enabled using the flag --preferTVarsWithStartValue=true. It solves this particular issue, but the model is still not working perfectly fine. Signal opAmp2.v_out shows some unexpected artefacts which need to be further investigated.

I will enable the option to privilege variables having start values, so that we can see its impact on the coverage test.

[Francesco Casella]

OK. Would you mind triggering the Jenkins job right away when you have changed the default, so we can see the impact of this single change?

[Lennart Ochel]

Francesco, please also note my comment on ​OMCompiler#3079.

[Francesco Casella]

Lennart, see my ​reply.

[Francesco Casella]

Lennart, as far as I understand, your ​PR #3079 pull request was good, as it fixed a number of problems in the testsuite. The remaining failures were probably caused by some bugs in the tearing algorithm, that should be addressed.

Overall, I would prefer to trade the many positive regressions on actual testsuite models, agains the negative regression on an academic test case. We should fix that, but we can do it later on.

Unfortunately, your PR ended up on a dead track because the OMCompiler project was archived.

Can you please resurrect it on the OpenModelica project, and push it in? I think we should have it in 1.14.0.

[Francesco Casella]

See ​https://github.com/OpenModelica/OpenModelica/pull/256

[Francesco Casella]

I checked the ​Jenkins report of the newInst testsuite after this build was merged. There are 29 improvements and 40 regressions. Unfortunately these overall figures are totally unreliable, because there are many errors caused by borderline timeouts, that cause random changes every time the testsuite is run. That's a quite annoying issue, see also #5184, but never mind.

Looking specifically into the most relevant libraries, we have

• BuildSysPro, Buildings, and IBPSA: (small) positive impact
• ModelicaTest 3.2.3: one more model passing verification
• ModelicaTest_trunk_cpp: five more models passing verification, two models going from Simulate to Compile (but they were not passing verification anyway)
• Modelica 3.2.3: two models going from Simulate to Verify, one going from Verify to Simulate, three going from Verify to Compile, one going from Simulate to Compile (this did not pass verification for an issue independent of the choice of tearing variables)
• ThermoPower: three models going from Simulate (there is no verification here) to Compile

Overall, we have mixed results: we solved the issue of the original ticket, we also got some more other models to work, but we also broke some models that were previously working.

The question is: is this inherently due to preferring variables with explicitly set start values for tearing, or is there still some bug there in the implementation that causes the new failures, e.g., because it doesn't recognized the variables with modified attribute correctly?

At this moment I couldn't say, and unfortunately I don't have time to look into this until next week.

@lochel, to help me with the analysis it would be very helpful if you could take the 4 MSL 3.2.3 examples that went from Verify to Compile (i.e., the simulation fails), and for each one get the list of tearing variables of the failing NLS, with their corresponding start values. Then, get the same list for the same systems with the previous version of OMC, which led to successful simulation, and post the two lists here.

By looking at them, we should probably get a good understanding of what happens and make the right decision.

Thank you!

[Lennart Ochel]

I guess it is mostly related to the optimization module wrapFunctionCalls. I had to deactivate it, because it caused problems combined with the tearing changes. This is something that needs to be investigated further.

[Francesco Casella]

Replying to lochel:

I guess it is mostly related to the optimization module wrapFunctionCalls. I had to deactivate it, because it caused problems combined with the tearing changes.

Aha. As I understand, wrapFunctionCalls is essential, also for performance reasons.

What kind of problems does this cause? We definitely can't ship 1.14.0 without wrapFunctionCalls.

[Francesco Casella]

Replying to lochel:

I guess it is mostly related to the optimization module wrapFunctionCalls. I had to deactivate it, because it caused problems combined with the tearing changes.

Do you mean removeEqualFunctionCalls? I can't find any mention of wrapFunctionCalls in the ​commit changes.

BTW, we already had some discussion about removeEqualFunctionCalls being obsolete in #4881 and #4952. There was no conclusion, maybe it's really time to sort that out.

[Vitalij Ruge]

Replying to lochel:

I guess it is mostly related to the optimization module wrapFunctionCalls. I had to deactivate it, because it caused problems combined with the tearing changes. This is something that needs to be investigated further.

You mean the reason why e.g. testAlgLoop5.mos failed, which was disabled.

My anaylsis for testAlgLoop5.mos:
For the analysis I use the simulation appi ( otherwise somebody can say it's because optimization is "implemented in a horrible way")

setCommandLineOptions("-d=dumpLoops");
getErrorString();

simulate(testAlgLoop5, tolerance = 1e-12);
getErrorString();

After ​https://github.com/OpenModelica/OpenModelica/pull/256 the model stop simulate.

The test work simulate with

omc testAlgLoop5.mos  --preOptModules+=removeEqualFunctionCalls

The output without removeEqualFunctionCalls

internal vars                                                                           
1: $cse3:VARIABLE(protected = true )  type: Real  unreplaceable
2: t4:VARIABLE(min = 1.0 max = 5.0 start = 5.0 )  type: Real                                               
3: $cse2:VARIABLE(protected = true )  type: Real  unreplaceable                         
                                                             
residual vars                                                                                           
1: t3:VARIABLE(start = 5.0 )  type: Real                                                       
2: t2:VARIABLE(min = -2.0 max = 0.5 start = 1.0 )  type: Real      
3: t1:VARIABLE(min = -0.4 start = 1.0 )  type: Real
                               
internal equation             
1/1 (1): $cse3 = exp(t1 + t2 + t3)   [unknown |0|0|0|0|]
2/2 (1): t3 = t1 + t2 + y + t4 + u3   [dynamic |0|0|0|0|]
3/3 (1): $cse2 = log(t4)   [unknown |0|0|0|0|]
                                 
residual equations         
1/1 (1): 0.0 = log(t1 + t2 + t3 + u + u2)   [unknown |0|0|0|0|]
2/2 (1): u1 = log((2.0 * t3 + 3.0 * t1 + t2) ^ 2.0)   [unknown |0|0|0|0|]                                               
3/3 (1): $cse3 = $cse2   [dynamic |0|0|0|0|]  

and with removeEqualFunctionCalls

internal vars                                                                                                                                                                                                      
1: t3:VARIABLE(start = 5.0 )  type: Real                                                                                                                                                                           
2: $cse3:VARIABLE(protected = true )  type: Real  unreplaceable                                                                                                                                                    
3: t4:VARIABLE(min = 1.0 max = 5.0 start = 5.0 )  type: Real                                                                                                                                                       
4: $cse2:VARIABLE(protected = true )  type: Real  unreplaceable                                                                                                                                                    
                                                                                                                                                                                                                   
residual vars                                                                                                                                                                                                      
1: t2:VARIABLE(min = -2.0 max = 0.5 start = 1.0 )  type: Real                                                                                                                                                      
2: t1:VARIABLE(min = -0.4 start = 1.0 )  type: Real                                                                                                                                                                
                                                                                                                                                                                                                   
internal equation                                                                                                                                                                                                  
1/1 (1): con = 2.0 * t3 + 3.0 * t1 + t2   [binding |0|0|0|0|]                                                                                                                                                      
2/2 (1): $cse3 = exp(t1 + t2 + t3)   [unknown |0|0|0|0|]                                                                                                                                                           
3/3 (1): t3 = t1 + t2 + y + t4 + u3   [dynamic |0|0|0|0|]                                                                                                                                                          
4/4 (1): $cse2 = log(t4)   [unknown |0|0|0|0|]                                                                                                                                                                     
                                                                                                                                                                                                                   
residual equations                                                                                                                                                                                                 
1/1 (1): 0.0 = log(t1 + t2 + t3 + u + u2)   [unknown |0|0|0|0|]                                                                                                                                                    
2/2 (1): $cse3 = $cse2   [dynamic |0|0|0|0|] 

At the end the size of the residual equation are change.

[Francesco Casella]

@vitalij, I don't really understand what removeEqualFunctionCalls does. I report here the code of ​testAlgLoop5:

model testAlgLoop5
  parameter Real one = 1;
  Real t1(start=1, min = -0.4) annotation(tearingSelect = always);
  Real t2(start=1,min=-2,max=1/2) annotation(tearingSelect = always);
  Real t3(start=5);
  Real t4(start=5, min=1, max=5);
  Real x(start=1,fixed=true,min=-3,max=3.0);
  Real y(start =1,fixed=true);
  Real z(start =-1,fixed=true, max=0.5, min = -1);
  input Real u(min=-5,max=5, start=1);
  input Real u1(min=-1,max=2, start=1);
  input Real u2(min=-1,max=2, start=1);
  input Real u3(min=-1,max=1, start=1);
  Real cost annotation(isLagrange = true);
  Real costM = sin(u*u1) annotation(isMayer = true);
  Real dummyCon = one*con annotation(tearingSelect = always);
  Real con(min=1, start=1, max = 2.5) = 2*t3 + 3*t1 + t2 annotation(isConstraint=true);
  Real conDer(min=-2, start=1, max = 2) = der(x) annotation(isConstraint=true);
  Real fcon(min=0,max=0) = 10*der(x) annotation(isFinalConstraint = true);
  Real fcon2(min=0) = u1 + u annotation(isFinalConstraint = true);
  Real fcon3(min=10,max=10) =10*( t1 + t2 + t3 + u3) annotation(isFinalConstraint = true);
  Real fcon4(min=10/2) = 10*u annotation(isFinalConstraint = true);
  parameter Real tgrid[:] = {i*50/250 for i in 1:250}  annotation(isTimeGrid = true);
equation
  t3 = t1 + t2 +y + t4 + u3;
  log(t1 + t2 + t3 + u + u2) = 0;
  log((2*t3 + 3*t1 + t2)^2) = u1;
  log(t4) = exp(t1 + t2 + t3);
  der(x) = x + t1*t2 + u*t3;
  der(y) = u2*x + y*z;
  der(z) = y;
  cost = cos(x)-sin(u*u1) + (u-1)^2 + u1^2 + (der(y)-sin(time))^2;
end testAlgLoop5;

I understand there is a strong component with unknowns t1, t2, t3, t4. When removeEqualFunctionCalls is deactivated, this loop is resolved by selecting t1, t2, t3 as tearing variables; addtionally, there are to $cse expressions, I guess due to wrapFunctionCalls. The first internal equation is solved for $cse3, the second for t4, the third for $cse2, and then there are three implicit equations to solve for t1, t2, t3.

When removeEqualFunctionCalls is turned on, t1, t2 only are selected as tearing variables. This means that the first internal equation has to be solved for t3 and the third for t4. But how can the first internal equation be solved for t3 if con is unknown? From what I understand of the structure of this problem, con only appears in the binding equation, so it must be solved once t1, t2, t3 are known. How can it be used in the torn equation before it has been computed?

[Vitalij Ruge]

removeEqualFunctionCalls replace

Real con(min=1, start=1, max = 2.5) = 2*t3 + 3*t1 + t2;

in

log((2*t3 + 3*t1 + t2)^2) = u1;

Following

log(con^2) = u1;

where con can be solved.

[Francesco Casella]

Thanks @vitalij for the explanation.

I've been trying for a long time to understand what removeEqualFunctionCalls does (see also #4881, #4952), maybe I'm getting closer to that. I beg your pardon if I cannot Read The Fine Source Code, but if you are not an OMC developer, understanding backend code written in MetaModelica is not trivial, unless someone guides you into it.

The ​documentation of this feature says:

removeEqualFunctionCalls (Detects equal function calls of the form a=f(b) and c=f(b) and substitutes them to get speed up.)

In this case we have

x = <expr>
log(<expr>²) = u1

being transformed into

x = <expr>
log(x²) = u1

which is clearly out of the scope of what the documentation says. If I try to reverse-engineer this behaviour, I'd say that what is actually done by this feature can be described as follows

when an equation x = <expr> exists, symbolically replace all other occurences of <expr> with x

Of course what is explained in the documentation is a special case, whereby <expr> is actually a function call.

Is this what removeEqualFunctionCalls does? If that is the case, there are two consequences

1. I guess we should rename it to removeEqualExpressions and update the documentation accordingly.

2. this functionality is not superseded by wrapFunctionCalls, which would only handle the special case of <expr> being a function call. This would sort out the questions raised in #4881 and #4952.

@lochel, @vitalij, can you please shed some light here, also regarding the mutual relationship between these two backend functions?

I hope this then helps understanding why removeEqualExpressions causes issue with the updated tearing algorithm, and allow to re-enable it. As I understand so far, both are needed, and I am not surprised that switching removeEqualExpressions off broke some models in the testsuite.

[Vitalij Ruge]

Some other effect is:
(Probably already known, but I wanted to emphasize it. If x replaced then x.start and x.nominal has impact of the loop. Special I think on Real der_x(nominal=1e-3) = der(x);.

[Lennart Ochel]

Thanks Francesco for the clarification; I agree completely. Actually, I haven't been aware of this undocumented functionality, and hence I thought the module is superseded.

I think it would be good to merge the functionality with one of the common sub-expression modules, and state the functionality clearly in the documentation. I propose to discuss it with the backend people in the next dev meeting.

[Francesco Casella]

OK, very good. So, maybe it's better to wait for the analysis of the recently broken models until this issue has been resolved - maybe they will heal automatically.

[Francesco Casella]

Summary of today's discussion: Lennart thinks that the additional feature brought in by removeEqualFunctionCalls is to perform CSE with expressions which are involved in function call arguments (not in general). If this is the case, then this could be integrated in the functionaly of comSubExp. He plans to do that as soon as he gets a working PC back.

[Francesco Casella]

@lochel, I guess it's more than that, and it's really more like I guessed in comment:14. Try the following test model

model M
  Real x, y, w, z, a;
equation
  a = time;
  a = x^2 + y^2 + w^2 - 3;
  x + 3*y - w  + z = 1;
  x - 3*y +2 -z = 2;
  x^2 + y^2 + w^2 -3 = z;
  annotation(__OpenModelica_commandLineOptions = "--preOptModules+=removeEqualFunctionCalls");
end M;

In this case, removeEqualFunctionCalls is actually able to single out the expression a = x^2 + y^2 + w^2 - 3 and use it to figure out that z = a, so it makes z an alias of a and that's it.

If you turn removeEqualFunctionCalls out this doesn't happen and z is computed as the torn variable of an implicit system. So, I think we agree removeEqualFunctionCalls does some CSE, but it's not limited to expressions appearing in function calls, which are absent in this example.

That said, I think moving it to comSubExp is fine, but we should understand exactly the scope and extent of what it does, not just make assumptions.

[Francesco Casella]

Ping :)

Is there any unanticipated difficulty in carrying this out?

[Francesco Casella]

removeEqualFunctionCalls was reactivated in ​PR #341. Thanks @lochel for sorting this out before your vacations :)

The ​impact on the libraries is apparently mixed (21 improvements, 22 regressions). However, I suspect some of the regressions are spurious, most of them impact on the cpp runtime, which is often a bit more fragile to this kind of changes (no idea why), but anyway, MSL 3.2.3 has three improvements and no regressions, which is good for our goal with 1.14.0. ThermoPower also has two more models getting back to simulate.

I think we can close this specific ticket because the original issue has been solved. I opened #5600 on the topic of removeEqualFunctionCalls, with target to 2.0.0.