OpenModelica inconsistent initial condistions for DAEs

I'm trying to use openmodelica to solve simple DAEs like for a planar pendulum. I am running into a simple problem that I don't now how to solve. I am imposing an initial condition for the x and y coordinate of the pendulum. x0 = 0.5 and y0 = 0 works since the length of the pendulum is 0.5. However, when I want to use x0 = -0.5, it breaks and tells me that the init conditions are incosistent but I don't understand why since x^2 + y^2 = L^2.

Error with negative x initial condition

I have also tried other negative init conditions like x0 = -0.013437983982246 and y0 = -0.499819387965786 but with the same error code. However when I make the x0 value positive (by removing minus-sign), modelica runs smoothly without any errors.

As you may have already found out, your model has a differentiation index higher than one. Unfortunately, solution algorithms like dassl and ida are only capable of solving index one problems. This is no shortcoming of dassl and ida but is a structural problem of high index.

How to find out the index of a DAE system? The index is the minimum number of times the algebraic equations of your DAE-system must
be differentiated with respect to time in order to receive a pure ODE system. In your case that has to be done three times. In order to check if your system is of index one, there is a simple algorithm: You can try to locate in each equation of your system exactly one "variable" that should be virtually calculated with that equation. "Variables" are all algebraic variables and the first order derivatives of the state variables (the states itself are no variables since they are given as initial values or provided by the integrator). If this works for a DAE system, then the index is one. Just repeat the differentiation step until you reach an index one problem.

For the numerical solution of high index problems several methods are available. Modelica introduces dummy variables (--> your favorite search engine), you can use projection methods and you can also transform your problem to a lower index by differentiation of the algebraic variables with respect to time.

The automatic introduction of dummy variables by openmodelica changes your equation system in the background, these changes are not visible (as long as you do not take a look at the generated c-code) and can and will have influence on the convergence of your nonlinear equation system. Just by chance I found out that after introducing the new variable Real L_test and the equation L_test = sqrt(x^2+y^2) your system solves for x0 = -0.5. Although this has nothing to do with the original equation, the additional equation changes Newtons's iteration path. It seems, that your error message is caused more or less because openmodelica did not converge due to the new, extended structure of the nonlinear equation system. Why did I introduce this new variable? See below...

How can we transform your problem to a lower index? There is only one algebraic equation

x^2 + y^2 = L^2;

which can be differentiated one time for receiving an index two problem

2*x*der(x) + 2*y*der(y) = 0;

and a second time for an index one problem:

2*der(x)^2 + 2*x*der(vx) + 2*der(y)^2 + 2*y*der(vy) = 0

(Of course the factor 2, for the sake of completeness, is not necessary in both equations)

If you use the first derivative (instead of x^2+y^2=L^2) your system is of index 2, Openmmodelica introduces one dummy variable and has no problem. You can also use the second derivative for index one without additional dummy variables. In both cases the system solves for your initial conditions x0 = -0.5 and x0 = 0.5.

Be aware of two points:

• After differentiating x^2 + y^2 = L^2 the physical meaning of this equation that the length is L is lost! You have to choose your initial conditions such that length-constraint is fulfilled.
• For the original DAE system the constraint on the length is always fulfilled. For the transformed system (index two or index one) this is not the case! How well the constraint is fulfilled depends on the integration accuracy. The longer the simulation time, the larger the error will be. That was the reason for introducing L_test = x^2 + y^2.

BTW: Your system has two degrees of freedom (x and y) and one constraint on the length L such that only one degree of freedom remains (which causes all the trouble). If you directly understand your system as having only one degree of freedom and expressing this by choosing polar coordinates, things are much easier.

PS: Next time please post the original model as source code, not as a picture.