Numerical overflow due to replacement of 1/tanh with cosh/sinh

[rfranke]

Some load models of the PowerSystems library adapt the load to make it independent of voltage. For instance PowerSystems.AC3ph.Loads.PQindLoad states:

  der(Z) = (pq/(pq*pq)*v2 * tanh(imax)/tanh(imax/V2_nom*v2) - Z) / tcst;

Only Z and v2 are time varying, whereas all other variables are parameters. OpenModelica replaces this with:

  der(Z) = (pq/(pq*pq)*v2 * tanh(imax)*cosh(imax/V2_nom*v2)/sinh(imax/V2_nom*v2) - Z) / tcst;

During the solution of the example PowerSystems.Examples.AC3ph.Inverters.InverterToLoad with DASSL it happens that v2 takes large values leading to an overflow of sinh. See the respective error message (sinh(...) becomes $cse11):

assert            | debug   | division leads to inf or nan at time 0.00374999, (a=inf) / (b=inf), where divisor b is: (pqLoad.pq[1] ^ 2.0 + pqLoad.pq[2] ^ 2.0) * $cse11

The overflow would not happen with tanh stated in the model as it is limited by 1. See e.g. https://en.wikipedia.org/wiki/Hyperbolic_functions:

[perost]

The code responsible for this is (note that the comment on the first line is wrong compared with what the code actually does):

OpenModelica/OMCompiler/Compiler/FrontEnd/ExpressionSimplify.mo

Lines 4786 to 4792 in 28cd7de

	// e1/tanh(e2) => e1*cos(e2)/sin(e2)
	case(_,op2 as DAE.DIV(ty),e1,DAE.CALL(path=Absyn.IDENT("tanh"),expLst={e2}),_,_)
	equation
	e3 = Expression.makePureBuiltinCall("sinh",{e2},ty);
	e4 = Expression.makePureBuiltinCall("cosh",{e2},ty);
	e = DAE.BINARY(e4,op2,e3);
	then DAE.BINARY(e1,DAE.MUL(ty), e);

I'm not sure what the purpose of this simplification rule is since it makes the expression bigger. Perhaps it enables further simplifications in some cases, but not in this case.

[rfranke]

This appears to go back to work by @vruge in 2014 first introducing e1/tanh(e2) => e1*coth(e2)
(OpenModelica/OMCompiler@aa6621c#diff-ba1b057b0ef11a721e80efc8709c6292974534491031177e8146309998c52648)
and then changing to e1/tanh(e2) => e1*cosh(e2)/sinh(e2)
(OpenModelica/OMCompiler@269a95a#diff-ba1b057b0ef11a721e80efc8709c6292974534491031177e8146309998c52648).

I think that this rule should be removed.

[rfranke]

There are more strange simplifications in ExpressionSimplify.mo:

  cos(e2)/tan(e2) => sin(e2)
  cosh(e2)/tanh(e2) => sinh(e2)

besides

  e1/tan(e2) => e1*cos(e2)/sin(e2)

[sjoelund]

The first two make a lot of sense though. You remove divisions by zero at least

[rfranke]

Commit a2907a1 fixes the first two.

[casella]

I completely agree that the simplification rule tanh(x) -> sinh(x)/cosh(x) is not a good idea, because of overflow. sinh(x) and cosh(x) grow exponentialy, while tanh(x) doesn't. Furthermore, tanh(x) is often used as a way to approximate step(x) with a smooth, analytic function, so it is quite likely that its argument becomes quite large.

I would definitely suggest to remove that rule. @rfranke, can you take care of that?

[rfranke]

OK, I merged PR #8006.

[casella]

PowerSystems.Examples.AC3ph.Inverters.InverterToLoad now runs successfully both with cpp runtime and c runtime, see regression report