Support the homotopy() operator for the generation of symbolic jacobians for NLS

[Francesco Casella]

Please run the MWE Test1 in the attached package with the declarative debugger. Among the initial equations, the analytic Jacobian is generated for the first nonlinear system, but not for the second, due to the use of the homotopy operator.

Setting --replaceHomotopy=actual, as in Test2, does not change the situation. This is strange, as I would have expected this flag to actually replace all instances of homotopy() with their actual expression - isn't that the case?

Setting -d=forceNLSanalyticJacobian, as in Test3, leads to the generation of the analytic jacobian for the second system:

(residual) homotopy(sin(y2), y2) - y2 + homotopy(sin(y1), y1) + y1 = 0
(residual) -1.0 + y1 + homotopy(sin(y1), y1) + y2 + homotopy(sin(y2), y2) = 0
(jacobian) $cse3 := cos(y1)
(jacobian) $cse4 := cos(y2)
(jacobian) $res_NLSJac453_1.$pDERNLSJac453.dummyVarNLSJac453 := y1.SeedNLSJac453 + $cse3 * y1.SeedNLSJac453 + y2.SeedNLSJac453 * (1.0 + $cse4)
(jacobian) $res_NLSJac453_2.$pDERNLSJac453.dummyVarNLSJac453 := y1.SeedNLSJac453 + $cse3 * y1.SeedNLSJac453 + $cse4 * y2.SeedNLSJac453 - y2.SeedNLSJac453

however, if I am not mistaken, while the residuals correctly contain the homotopy() operator, the jacobian seem to be computed assuming lambda = 1, and is thus wrong in general.

Differentiating the homotopy operator is no big deal, as the homotopy an differentiation operators are commutative. Is it possible to implement that right away?

[Andreas Heuermann]

I'm looking into this and see what we can do.

Let's start with some more analysis. Btw. @casella I like your MWE.

---

-d=symjacdump for Test1

Simplified equation system is looking good.

Create symbolic Jacobians from:
Independent Variables
========================================
1: y1:VARIABLE()  type: Real
2: y2:VARIABLE()  type: Real
Dependent Variables
========================================
1: $res_LSJac3_1:VARIABLE()  type: Real
2: $res_LSJac3_2:VARIABLE()  type: Real
Basic equation system:
differentiated equations
========================================
1/1 (1): $res_LSJac3_2 := 2.0 * y1
2/2 (1): $res_LSJac3_1 := -1.0 + 2.0 * (y1 + y2)
related variables
========================================
1: $res_LSJac3_2:VARIABLE()  type: Real
2: $res_LSJac3_1:VARIABLE()  type: Real
known variables
========================================
Symbolic Jacobian:

unknown partition
========================================

Variables (2)
========================================
1: $res_LSJac3_2.$pDERLSJac3.dummyVarLSJac3:VARIABLE()  type: Real unreplaceable
2: $res_LSJac3_1.$pDERLSJac3.dummyVarLSJac3:VARIABLE()  type: Real unreplaceable


Equations (2, 2)
========================================
1/1 (1): $res_LSJac3_2.$pDERLSJac3.dummyVarLSJac3 = 2.0 * y1.SeedLSJac3   [dynamic |0|0|0|0|]
2/2 (1): $res_LSJac3_1.$pDERLSJac3.dummyVarLSJac3 = 2.0 * (y1.SeedLSJac3 + y2.SeedLSJac3)   [dynamic |0|0|0|0|]

And for the actual system we have:

Create symbolic Jacobians from:
Independent Variables
========================================
1: y1:VARIABLE()  type: Real
2: $cse5:VARIABLE(protected = true )  type: Real unreplaceable
3: $cse6:VARIABLE(protected = true )  type: Real unreplaceable
4: y2:VARIABLE()  type: Real
Dependent Variables
========================================
1: $res_NLSJac5_1:VARIABLE()  type: Real
2: $res_NLSJac5_2:VARIABLE()  type: Real
3: $res_NLSJac5_3:VARIABLE()  type: Real
4: $res_NLSJac5_4:VARIABLE()  type: Real
Basic equation system:
differentiated equations
========================================
1/1 (1): $res_NLSJac5_4 := -1.0 + y1 + $cse5 + y2 + $cse6
2/2 (1): $res_NLSJac5_3 := sin(y1) - $cse5
3/3 (1): $res_NLSJac5_2 := y1 + $cse5 + $cse6 - y2
4/4 (1): $res_NLSJac5_1 := sin(y2) - $cse6
related variables
========================================
1: $res_NLSJac5_4:VARIABLE()  type: Real
2: $res_NLSJac5_3:VARIABLE()  type: Real
3: $res_NLSJac5_2:VARIABLE()  type: Real
4: $res_NLSJac5_1:VARIABLE()  type: Real
known variables
========================================
Symbolic Jacobian:

unknown partition
========================================

Variables (4)
========================================
1: $res_NLSJac5_4.$pDERNLSJac5.dummyVarNLSJac5:VARIABLE()  type: Real unreplaceable
2: $res_NLSJac5_3.$pDERNLSJac5.dummyVarNLSJac5:VARIABLE()  type: Real unreplaceable
3: $res_NLSJac5_2.$pDERNLSJac5.dummyVarNLSJac5:VARIABLE()  type: Real unreplaceable
4: $res_NLSJac5_1.$pDERNLSJac5.dummyVarNLSJac5:VARIABLE()  type: Real unreplaceable


Equations (4, 4)
========================================
1/1 (1): $res_NLSJac5_4.$pDERNLSJac5.dummyVarNLSJac5 = y1.SeedNLSJac5 + $cse5.SeedNLSJac5 + y2.SeedNLSJac5 + $cse6.SeedNLSJac5   [dynamic |0|0|0|0|]
2/2 (1): $res_NLSJac5_3.$pDERNLSJac5.dummyVarNLSJac5 = cos(y1) * y1.SeedNLSJac5 - $cse5.SeedNLSJac5   [unknown |0|0|0|0|]
3/3 (1): $res_NLSJac5_2.$pDERNLSJac5.dummyVarNLSJac5 = y1.SeedNLSJac5 + $cse5.SeedNLSJac5 + $cse6.SeedNLSJac5 - y2.SeedNLSJac5   [dynamic |0|0|0|0|]
4/4 (1): $res_NLSJac5_1.$pDERNLSJac5.dummyVarNLSJac5 = cos(y2) * y2.SeedNLSJac5 - $cse6.SeedNLSJac5   [unknown |0|0|0|0|]

But we don't have a system for the equations with the homotopy operator.

[Andreas Heuermann]

For Test2 there is not a single change in the Jacobians.

Apparently we have two ways to replace the homotopy operator, ​elabBuiltinHomotopy in the old Frontend and ​replaceHomotopyWithSimplified in the BackEnd.

I checked which function to replace the homotopy was called and only the Backend function is called. But it is not checking the value of the flag replaceHomotopy and always replacing with the simplified system.

Before:

Equations (2, 2)
========================================
1/1 (1): y1 + homotopy(sin(y1), y1) + homotopy(sin(y2), y2) - y2 = 0.0   [dynamic |0|0|0|0|]
2/2 (1): y1 + homotopy(sin(y1), y1) + y2 + homotopy(sin(y2), y2) = 1.0   [dynamic |0|0|0|0|]

After:

Equations (2, 2)
========================================
1/1 (1): y1 + y1 + y2 - y2 = 0.0   [dynamic |0|0|0|0|]
2/2 (1): y1 + y1 + y2 + y2 = 1.0   [dynamic |0|0|0|0|]

EDIT:
replaceHomotopy is completely ignored in the newFrontend. Of course we always need replaceHomotopyWithSimplified for the initialization with lambda=0.

[Andreas Heuermann]

@perost Can we add something similar to ​elabBuiltinHomotopy to the new frontend? It sould be very easy, but I don't know what the best location would be.

[Per Östlund]

Replying to AnHeuermann:

@perost Can we add something similar to ​elabBuiltinHomotopy to the new frontend? It sould be very easy, but I don't know what the best location would be.

Sure, shouldn't be much of a problem. I'll see if I can add it.

[Per Östlund]

I made a PR which hopefully fixes the frontend issue, see ​#7069.

It's not quite equivalent to the old frontend because it types the whole call first and then simplifies it, while the old frontend does the reverse. I think it's better to make sure the call is actually valid even if it's going to be simplified, and it also makes the code simpler.

[Andreas Heuermann]

I added a rule for differentiating the homotopy operator in ​PR#7070.

I suspect I need to update some tests and maybe the runtime as well, but the differentation module for homotopy is now looking good.

[Francesco Casella]

Thank you folks, I really appreciate the quick fix.

I'll wait for the two PRs to be merged in and then I'll try them out on the original problem, which is way more involved than the MWE. @AnHeuermann, I am happy you liked it, William of Ockham would be proud of how I apply his famous razor to OMC debugging :)

[Francesco Casella]

@AnHeuermann, based on the comments in ​PR 7070 I updated the attached MWE package to actually use the sparse Kinsol/KLU solver, as will be the case for the actual application (large-scale power grid modelling) that motivated this issue.

More specifically, Test6 directly forces to use Kinsol as nonlinear solver and to use KLU as solver for Kinsol. Please make sure that all the required code is generated to support that case.

BTW, I thought I would get the same result by resetting the min size and max density, see Test5. For some reason, instead, that test case still uses the default dense solver. Is that a but, or do I miss something?

[Andreas Heuermann]

You still need -d=forceNLSanalyticJacobian for a symbolic Jacobian to be generated, but it should work.
I was looking at old code when I couldn't find the Jacobian in the generated code.

I corrected the failing tests and added a test as well. It's a bit to verbose for my taste, but it has to do for now.

[Francesco Casella]

Replying to AnHeuermann:

You still need -d=forceNLSanalyticJacobian for a symbolic Jacobian to be generated, but it should work.

Why should I need that? The explanation of the flag states:

Forces calculation (of) analytical jacobian also for non-linear strong components with user-defined functions.

Now, homotopy() is a built-in function, not a user-defined one. It's the same as sin() and cos(). As I understand, the rationale of that flag is to avoid the use of AD to generate the "differentiated algorithms" in user-defined functions, maybe because it's not too reliable or potentially inefficient in some cases. There's no such need with homotopy(), which is trivially differentiated.

From an end-user perspective, there's nothing fancy in the PowerGrids models, except the need of using daeMode. They should run efficiently out of the box without the need of magic incantations in the form of arcane debug flags.

@AnHeuermann, could you please also fix that? It's not urgent, but it would be good to have in 1.17.0, which people will use for a while to run PowerGrids models.

In fact, even the need of using daeMode should be detected automatically, since the DAE has a sparse Jacobian, while the ODE has a dense one. I opened #6340 on this topic.

[Andreas Heuermann]

I'm not sure what the criteria are to generate a symbolic Jacobian for non linear systems.
@kabdelhak Do you know that?

[Francesco Casella]

I remember discussing this with Willi Braun before he left. What I remember is what I reported in comment:13, i.e., symbolic Jacobians are generated by default in case only built-in functions show up in the nonlinear equation system, while they are not if other functions show up, unless -d=forceNLSAnalyticJacobian is set.

I checked the ​calculateStrongComponentJacobians module and I found the ​checkForSymbolicJacobian function, which I guess is where the choice between symbolic and numerical is made. I'm not really able to fully follow the logic of this MetaModelica code, but I guess the ​hasEqnNonDiffParts function is the culprit. Shouldn't we just remove ​this line?

[Andreas Heuermann]

That was indeed the solution. Thanks!
I'll update my PR.

[Francesco Casella]

Wow, I lost the opportunity of committing my first line of code to the codebase, too bad :P

I hope you can manage to merge next week at the latest, so I can get it in the nightly build.

Thanks!

[Francesco Casella]

Problem solved by ​PR 7070.