How to add “nonnegative” to DAE

I have solved the DAE solution using ODE15i, but the data in the result is negative. I need to solve for the angle and I need the angle and the velocity to be non-negative. But you can not apply the "non-negative" option when solving DAE problems. And I have tried to reduce "ReTol" and "AbsTol", but I still do not get ideal results. Do you have any other plans?

[Torsten]

Totally differentiate the algebraic equations twice with respect to t and set up your system as an ODE system of 16 ODEs instead of a DAE system.

[Linux Jia]

I have tried this method, but because the system is too complicated, the result will not satisfy the original algebraic equation.

[Torsten]

Which of the variables become negative and should come out to be positive ? All angles and angular velocities ? Or with which angles and angular velocities does it start that they become negative ? How big is "negative" ? -1e-6 or bigger ?

[Linux Jia]

The angular velocity cannot be negative, the initial value is +0.1, it becomes negative after 50 data points, and continues to become negative, and the angular velocity shows a linear downward trend. The physics system requires both angles and velocities to be non-negative.

[Torsten]

Do your initial values for phi1,...,phi8,phidot1,...,phidot8 satisfy your algebraic constraints (also in differentiated form) and your implicit differential equations ?

[Linux Jia]

Verify that the initial values of the angle and angular velocity satisfy the equations by decic.

[Torsten]

You mean you verified this using "decic" ?

[Linux Jia]

Yes.

[Torsten]

Just out of interest:

How is it possible to fix initial conditions for phi1,...,phi8,phi1dot,...,phi8dot that guarantee that angle and velocity remain positive throughout the integration and additionally satisfy the algebraic constraints ?

[Linux Jia]

Decic can ensure that the initial values of the angle and angular velocity satisfy the equation, but may be negative. A reasonable initial value can be guaranteed to satisfy the equation and may be non-negative after being verified by decic. And now my situation is like this, but it is not always satisfied that the velocity is non-negative during the integration.

[Torsten]

But you are sure that for your system, there should be initial values for phi and phidot that guarantee that angle and angular velocity remain positive during integration ? If yes: why can you be sure about this ? From the physics of the problem ?

[Linux Jia]

Because the physical system has a ratchet mechanism, the angle is only allowed to turn in the positive direction. So the angle cannot be negative, and the angular velocity cannot be negative.

[Torsten]

Ok, then the initial conditions are either not unique so that they don't represent the physical state correctly or there is a mistake in your code.

[Linux Jia]

The initial values are correct and there is no error in the equation code. There's a reason for that. It is because I have written a physical system of 6 angles by differentiating the algebraic equations, without errors, and the results obtained satisfy the algebraic equations without differentiation. And the current equation code has not changed, because my program is automatic, just need to set the value of "n" to describe the system of different magnitudes. This can prove that my code is correct, and the method of setting the initial value is correct.

[Torsten]

And what makes the results "unphysical" then in your opinion ?

[Linux Jia]

Maybe the function of matlab can't make the DAE system non-negative? Because conventional solvers such as ode45 can impose non-negative.

[Torsten]

In all cases when I used the non-negative option, this option couldn't help the solution to remain positive.

The reason is that if the equations don't permit something else but a negative solution, also this option cannot prevent the solution to become negative.

In my opinion, it's very unlikely that amplified rounding errors can lead to massive negativity. And you said you already reduced the default tolerances of the implicit integrator. That's all you can do from the solver's side. Or you could try a different implicit integrator, e.g. dassl, radau5 or similar ones.

[Jeffrey Clark]

Angles are funny things:

1. Generally every negative has its positive when dealing with trig functions, so are you also saying that combinations of angles that would result in different signs of the trig function combinations (or the one) is invalid?
2. Is it possible that your solution that finds negative angles and velocity isn't really valid when concidering the angle may be pointing directly opposite what you would expect, which would also therefore have a negative velocity and be a valid mathematical solution, but would need to be interpreted correctly to agree with your apperatus?
3. Is the system of equations and functions such that there has to be some 'unwrap'-like done?

[Linux Jia]

I think this is a multi-solution problem and hopes to get the results of the physical system. Alternatively, can the ODE solverator process the segmentation function?

[Torsten]

What is the segmentation function ?

[Linux Jia]

piecewise equation

[Torsten]

Is this a different question or does it have to do with the problem from above ? Since I don't see where piecewise function come into play here.

[Linux Jia]

That is, when the angular velocity is negative, assign its value to 0.

[Torsten]

You can try, but if the negative value is more than just a slight numerical error that is corrected in the next iteration, I don't think this will lead you somewhere.

[Linux Jia]

I'm not rambling, two workarounds have been suggested in matlab answers before. The first is to reduce the integration tolerance, and the second is to add judgment in the integrator. If a variable is negative, set its value to a very small value. Also, this person seems to be you. ^-^

https://ww2.mathworks.cn/matlabcentral/answers/161302-how-to-solve-dae-system-with-non-negativity-constraint?s_tid=srchtitle

[Bruno Luong]

It looks to me that the system might have multiple solutions, and it is hard to catch the solution that is positive.

I would suggest then to start with a positve state phi_{k=0), then slowly iterate until your target DAE meets. Here is a pseudo code:

phi_{k=0) >= 0; % and also the 3 last equations
alpha = some value in (0,1)
for k=1,2,...
    Qp = ComputeQ(phi{k-1}) 
    % Use phi{k-1} as first guess solve DAE with RHS set to
    d(xxx)/dt ... = (1-alpha)*Qp + alpha*Q;
    if phi_{k} < 0
        alpha = (0+alpha)/2; % reduce alpha;
    else
        alpha = (1+alpha)/2; % eventually increase alpha
    end
    if alpha ~= 1
        % we solve the DAE approximately
        break
    end
end 

Answer 1

[Linux Jia]

I added a judgment statement to the handle function of the solver, and if the angular velocity is negative, set its value to an extremely small value. The program can continue to run, but when the time of the program reaches a certain point, it prompts that the step size needs to be reduced to meet the integration tolerance requirements. So how should this situation be handled? If the last numerical solution of this integration is used as the guessed initial value of the next integration, it will be found that its initial value has changed after passing decic. As a result, the junction of the two integrals will not be smooth. So, is there any good way to deal with it?

[Torsten]

From the time instant when you artificially reset the angular velocity to a small value, the numerical values you receive for the solution variables will no longer be solutions of the DAE-system. That's why decic produces a change in the initial values when you try to restart the solver. I think you will have to accept this.

Answer 2

[Linux Jia]

The integral does not converge because the mathematical model is wrong.

[Torsten]

And how did you find it out ?