Solving equations with a division by zero

[anonymous]

I attached two examples which are current not simulating in OpenModelica. I know these examples should work as they are simulating in Dymola and they have been working in a previous version two years ago (when I developed similar models). Please load the system library ElectroMechanicalDrives to investigate the following two examples.

TestElectroMechanicalDrives

When compiling this example the following error message is generated:

stdout | unknown | <p>DASKR--  AT T (=R1) AND STEPSIZE H (=R2) THE                                    <br>
      In above,  R1 =   0.0000000000000E+00   R2 =   1.6518123698891E-13<br>
DASKR--  NONLINEAR SOLVER FAILED TO CONVERGE                                    <br>
DASKR--  REPEATEDLY OR WITH ABS(H)=HMIN                                         <br>

stdout | error | <p>Simulation process failed. Exited with code 248.</p>
stdout | info | <p>The initialization finished successfully without homotopy method.</p>
stdout | warning | <p>The corrector could not converge.</p>
stdout | warning | <p>can't continue. time = 0.000000</p>
stdout | info | <p>model terminate | Integrator failed. | Simulation terminated at time 0</p>

When I remove the inertia, the model simulates OK. This I do not understand...

TestElectroMechanicalDrives2

When compiling this example the following error message is generated:

stdout | error | <p>Simulation process failed. Exited with code 1.</p>
stdout | info | <p>The initialization finished successfully without homotopy method.</p>
assert | debug | <p>division by zero at time 0, (a=0) / (b=0), where divisor b expression is: inertia.w</p>
assert | debug | <p>division by zero at time 0, (a=0) / (b=0), where divisor b expression is: inertia.w</p>
assert | debug | <p>division by zero at time 0, (a=0) / (b=0), where divisor b expression is: inertia.w</p>

This issue seems to be related with speed. If I change the offset of the ramp to 1E-3 the example works fine.

Summary

The problem seems to be caused by ElectroMechanicalDrives.Components.Rotational.ConstantEfficiency. However, this example should be working and has been working some time ago. As a workaround I wonder if there exists an alternative formulation of the efficiency equations to work with OpenModelica.

[Christian Kral <dr.christian.kral@…>]

Example 1

[Christian Kral <dr.christian.kral@…>]

Example 2

[Christian Kral <dr.christian.kral@…>]

Reported by Christian Kral, referring to

OMEdit 1.13.0~dev-1-g933120e
Connected to OpenModelica 1.13.0~dev-25-g7b13b9e
Linux Mint 18.2

[Christian Kral <dr.christian.kral@…>]

OK, I have some new information here: I traced back the problem to the attached little package. Please consider the attached example Efficiency.Drive. In this example the efficiency model Efficiency.ConstantEfficiency is used.

CASE 1

In Efficiency.ConstantEfficiency the constant efficiency is considered by the following core equations:

  Modelica.SIunits.Power power_a(start=0) = flange_a.tau * w 
    "Power input of flange_a";
  Modelica.SIunits.Power power_b = flange_b.tau * w 
    "Power input of flange_b";
equation
  flange_b.tau = smooth(0, if power_a > 0 then -flange_a.tau * efficiency 
                                          else -flange_a.tau / efficiency);

In this implementation the attached Efficiency.Drive works OK.

CASE 2

But, if I change the equation in the efficiency model to (the commented line)

  power_b = smooth(1, if power_a > 0 then -power_a * efficiency 
                                     else -power_a / efficiency);

the example Efficiency.Drive is not working. The problem with CASE 2 seems to be, that for zero speed (w=0) torque is not calculated correctly since "only" a power equation is stated. In CASE 2 the power euqations

power_a = flange_a.tau * w
power_b = flange_b.tau * w

are currently not automatically divided by w in OpenModelica...

[Christian Kral <dr.christian.kral@…>]

Test example tracing the problem

[wbraun]

After all we end up here with the equation:

power_a = flange_a.tau * w

which is solved for flange_a.tau and this results into a division by zero since w=0 and we get the reported issue.

The expected solution is flange_a.tau=0 for time=0:0.1, right?

When we would handle this equation as a linear system then our solution would be zero.

The question is how detect such cases, is the occurrence of a division enough as criteria?

[dr.christian.kral@…]

Re-added the two examples with intertia (as intended). When using the actual release version 1.0.0 of ElectroMagneticDrives, both examples simulate fine.

TestElectroMechanicalDrives

TestElectroMechanicalDrives2

Yes the issue is that the equation

power_b = smooth(1, if power_a > 0 then -power_a * efficiency 
                                     else -power_a / efficiency);}}}

shall be divided by w if w=0, reducing the equation to

flange_b.tau = smooth(0, if power_a > 0 then -flange_a.tau * efficiency 
                                          else -flange_a.tau / efficiency);

as a consequence of

power_a = flange_a.tau * w
power_b = flange_b.tau * w

[dr.christian.kral@…]

New version of example 1 with inertia

[dr.christian.kral@…]

New version of example 2 with inertia

[dr.christian.kral@…]

Using the implementation of version 0.13.0 of ElectroMechanicalDrives

power_b = smooth(1, if power_a > 0 then -power_a * efficiency 
                                     else -power_a / efficiency);}}}

the examples to not work.