Opened 5 years ago

Last modified 5 years ago

#5624 new discussion

if else ODE and Algebraic Equation

Reported by: Matheus Campanini Mughrabi <matheus.mughrabi@…> Owned by: Lennart Ochel
Priority: high Milestone: Future
Component: Backend Version: v1.13.2
Keywords: If Else ODE Algebraic Cc: Karim Adbdelhak, Andreas Heuermann

Description

Hello, everyone!

I am facing a problem with a model I`m developing in OpenModelica and I created I simpler example to show the issue:

model teste

Real x;

equation

if time<1000 then

der(x) = 10;

else

x = -5;

end if;

end teste;

So what happens is that when "time = 1000" the simulation stops, I imagine it`s a discontinuity issue. Does anybody know how to solve this problem?

Thank you!

Change History (6)

comment:1 by Karim Adbdelhak, 5 years ago

I don't know what you are trying to achieve, but following code does what i think was your intention:

model teste
    Real x;
equation
    if time<1000 then
      der(x) = 10;
    else
      der(x) = 0;
    end if;
    
    when time>1000 then
      reinit(x, -5);
    end when;
end teste;

Try to make sure that each branch of an if equation can be solved for the same variables. If you want to change a state discontinuously use the reinit() operator together with when.

Last edited 5 years ago by Karim Adbdelhak (previous) (diff)

in reply to:  1 comment:2 by Karim Adbdelhak, 5 years ago

Replying to Karim.Abdelhak:

Try to make sure that each branch of an if equation can be solved for the same variables. If you want to change a state discontinuously use the reinit() operator together with when.

Actually it is not really necessary to have each branch be solvable for a the same variables, it can be simulated as a residual. It just results in better runtime code if you make sure that they do. But when trying to make a state constant at a given event, you have to set the derivative to zero and use reinit.

comment:3 by Francesco Casella, 5 years ago

Cc: Karim Adbdelhak Andreas Heuermann added
Type: defectdiscussion

It is common wisdom in the community that Modelica tools cannot handle systems with variable number of state variables, as will be the case here, see, e.g., Dirk Zimmer's PhD dissertation. However, I failed to find any such indication in the Specification, which instead mentions explicitly "variable structure" DAEs. So, it seems that the example model is perfectly legal, from the point of view of the Specification.

@Karim, @lochel, @AnHeuermann, do you have any comments on this?

in reply to:  3 comment:4 by Karim Adbdelhak, 5 years ago

Replying to casella:

It is common wisdom in the community that Modelica tools cannot handle systems with variable number of state variables, as will be the case here, see, e.g., Dirk Zimmer's PhD dissertation. However, I failed to find any such indication in the Specification, which instead mentions explicitly "variable structure" DAEs. So, it seems that the example model is perfectly legal, from the point of view of the Specification.

@Karim, @lochel, @AnHeuermann, do you have any comments on this?

A correct approach for the compiler would be to formulate the if expression implicitly (this is actually already done in the backend):

0.0 = if time < 1000 then der(x) - 10 else x + 5; 

now we can see that the second branch is structurally singular in the sense of index reduction and has to be differentiated (however this should be detected, seems quite complicated for more difficult if equations):

0.0 = if time < 1000 then der(x) - 10 else der(x); 

Now it is perfectly legal, but the original constraint equation of the second if branch has to be made a reinit in a when, since it only holds at a specific point in time:

   when time>1000 then
      reinit(x, -5);
    end when;

In theory this is doable, but i don't think it is worth the effort. Even though the specification does seem to allow models such as the proposed one, i would rather update the specification to not allow it.

By the way: Dymola also cannot handle such models.

Last edited 5 years ago by Karim Adbdelhak (previous) (diff)

comment:5 by Francesco Casella, 5 years ago

I opened a discussion on the Modelica Specification issue tracker, see #2411. Please follow up there.

comment:6 by anonymous, 5 years ago

I actually figured out a way to solve the problem by using a "high-pass filter" such as follows:
if time<1000 then
der(x) = 10;
else
der(x) = B*(-5 - x);
end if;

where B is a parameter that you fix, such as B=1000 (the higher it`s the less it changes the dynamics of the system).
This is actually a solution that is commonly used by programmers of several languages.

Note: See TracTickets for help on using tickets.