Simulation of Modelica.Electrical.Analog.Examples.OpAmps.Subtracter goes astray

[casella]

Please check the simulation of Modelica.Electrical.Analog.Examples.OpAmps.Subtracter.bmo: the vOut.v should just be the diffence vIn1.v-vIn2.v, which is a sinusoidal signal; instead, it veers off its course multiple times.

[jClugstor]

Solving this with FBDF I get perfect results, so I don't think it's a model issue. It might just be that the default solver for OrdinaryDiffEq is not very good on this problem.

[casella]

Traditionally, the models in the Modelica Standard Library have all been simulated with DASSL (a variable-order, variable step-size BDF solver), after transforming them into explicit ODEs. Maybe not the optimal choice for performance (some models run faster with other solvers, of course including explicit solvers for non-stif problems), but always reliable.

Is something similar available in Julia, that we can use as a default for the whole testsuite?

[jClugstor]

There is a DASSL implementation in the Differential equations suite, so we should be able to just use that.

@AayushSabharwal do you know if it's possible to use DASSL.jl on problems generated from MTK? It looks like it only works on DAEProblems from reading the documentation.

[AayushSabharwal]

MTK can generate DAEProblems just fine. I'm not sure if DASSL solves initialization, though

[casella]

In Modelica tools, DASSL is often used also for ODEs. I guess mostly for historical reasons, 25 years ago it was the only open-source robust stiff adaptive solver with event detection that was available. Still difficult to beat in terms of robustness.

Initialization problems are solved separately by ad-hoc generated code.

[casella]

I'm not sure about the details of FBDF, but it's probably not too different from DASSL. I guess the problem here is with the event-handling part. The OpAmp models have high gain and voltage saturation, this is known to have caused problems in the past.

[AnHeuermann]

While parsing the Base Modelica file we get warnings

Model:   Modelica.Electrical.Analog.Examples.OpAmps.Subtracter
┌ Warning: Did not converge after `maxiters = 100` substitutions. Either there \
│ is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics /home/runner/.julia/packages/Symbolics/eUt8E/src/variable.jl:443
┌ Warning: Did not converge after `maxiters = 100` substitutions. Either there \
│ is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics /home/runner/.julia/packages/Symbolics/eUt8E/src/variable.jl:443

See Modelica.Electrical.Analog.Examples.OpAmps.Subtracter_parsing.log.

What do these mean?

And during solving with Rodas5P we get

Model:   Modelica.Electrical.Analog.Examples.OpAmps.Subtracter
┌ Warning: Verbosity toggle: rosenbrock_no_differential_states 
│  Rosenbrock methods on equations without differential states do not bound the error on interpolations.
└ @ OrdinaryDiffEqCore /home/runner/.julia/packages/OrdinaryDiffEqCore/mNgXc/src/solve.jl:163

See Modelica.Electrical.Analog.Examples.OpAmps.Subtracter_sim.log.

This explains why the result is so bad. Forcing a smaller fixed step size improves the result:

I'll switch to FBDF in the test script to keep algebraic variables under an error bound. We could also try to generate a DAEProblem and use IDA at some point. But this is simpler to do right now.

[ChrisRackauckas]

Or use Rodas5Pr

[ChrisRackauckas]

I'm not sure about the details of FBDF, but it's probably not too different from DASSL. I guess the problem here is with the event-handling part. The OpAmp models have high gain and voltage saturation, this is known to have caused problems in the past.

Did you actually just run DASSL instead of speculating? It's like 2 lines of code so that would clear things up. Pkg.add("DASKR"), then just use daskr() as the solver for a DAEProblem. If you suspect it's the solver, just do those two lines and turn speculation into something provable. Most likely it's not that, it's a problem formulation thing though, but feel free to prove it wrong.

There is a DASSL implementation in the Differential equations suite, so we should be able to just use that.

There is also DASSL.jl too. It's a direct translation of DASSL, but DASKR is generally a better thing to do (it's just DASSL with the linear solver swap) and when Modelica folk say they use DASSL, generally they mean they are actually using DASKR (at least when I looked at the OpenModelica binaries and such, "DASSL" actually pointed to the DASKR Fortran binary, which is literally what's wrapped in DASKR.jl, so it's a direct 1-1 comparison while the DASSL.jl is a translation which could have another variable involved).

Traditionally, the models in the Modelica Standard Library have all been simulated with DASSL (a variable-order, variable step-size BDF solver), after transforming them into explicit ODEs. Maybe not the optimal choice for performance (some models run faster with other solvers, of course including explicit solvers for non-stif problems), but always reliable.

DASSL is generally not a great choice in this day and age because while it's a BDF method, it's not a fixed leading coefficient BDF method, so it loses stability in comparison to FLC implementations due to its quasi-constant constraints. That's why FBDF, DFBDF, and IDA are generally preferred. Dymola has a "modified DASSL", which is able to handle some significant things that DASSL on its own cannot. It's not public code so the exact details of what is modified are not available, but (a) it's these modifications that make it so DASSL in Dymola handles what Sundials does not (which can be tested using its FMU builder, since the FMU builder doesn't include the modified DASSL), and (b) from what I have tested modern FBDF handles all of those special cases now. So while I can only speculate as to what was modified, I believe FBDF has all of those modifications (and more) by now, which is why it ends up strictly more stable than IDA in many of these cases. OpenModelica does not use a modified DASSL (DASKR), so the common cases where DASSL fails but IDA succeeds applies in OMC as well, and OMC's default seems to be a bit of legacy voodoo in the community thinking DASSL is the best when in reality, DASSL is the best choice in Dymola simply because it's not using the normal DASSL and OMC probably shouldn't follow that by defaulting to normal DASSL.

But anyways, tl;dr FBDF on most recent versions has a bunch of deep tricks that should be matching the Dymola tricks while also being FLC, which gives best of both worlds. Then with a nonlinear solver swap you can also make it do some wild things if you really want.

[AnHeuermann]

I played around with DASKR.jl and it is a bit tricky to use. Unless I used it wrong there is no dense output, so I had to use linear interpolation. Also initializing the DAEProblem seems to be tricky. I used the ODEProblem for that.

For Modelica.Electrical.Analog.Examples.CharacteristicIdealDiodes both Rodas5Pr and DASKR get to a decent solution:

The question is what we want to showcase in BaseModelicaLibraryTesting.jl. I think the greatest advantage of going Modelica -> Base Modelica -> MTK is the easy integration of all the available tools from the Julia ecosystem. So going with FBDF sounds like a good demonstration of the usage of modern solvers.

From my point of view this issue can be closed. The generated MTK model is correct. We just need to simulate it correctly for the coverage report.

[ChrisRackauckas]

Also initializing the DAEProblem seems to be tricky.

That should be handled in the most recent version, if not it would be good to get an MWE to play with.

[casella]

@ChrisRackauckas I'm not sure why you closed this ticket, but the problem is still there.

I see several options to solve it:

• Fix Rodas5P
• Change the default solver to Rodas5Pr or DASKR in BaseModelica.jl
• Change the solver from the default in our CI

What would you recommend?