Improve zero crossing functions by using continuous expressions

[Francesco Casella]

Zero crossing functions are currently implemented with discontinuous outputs. Conceptually speaking, x > 2 is mapped into

zcf = noEvent(if x > 2 then 1 else -1);

Solvers try to find the event time using iterative algorithms. IDA uses a modified secant method. With this implementation, it reduces to the bisection method, which halves the error at each iteration, or equivalently gains 3 decimal digits every 10 iterations.

If we had a continuous zcf like

zcf = x - 2

the secant algorithm would converge in a single step.

The question of course is how to handle more complicated conditions such as

if x > 2 and y < 3 and isActive then ...

where x and y are Real variables and isActive is a boolean variables. A trivial discontinuous implementation is

zcf = noEvent(if x > 2 and y < 3 and isActive) then 1 else -1)

but it won't be very efficient. A better one could be to consider

zcf1 = x - 2;
zcf2 = y - 3;

separately, or possibly to consider their product

zcf = (x - 2)*(y - 3)

which will change the sign whenever one of the two conditions switches from false to true or vice-versa.

[Francesco Casella]

When this is done, also remember to update spatialDistribution and delay, see ticket:5281#comment:18

[Andreas Heuermann]

Delay doesn't have a zero crossing function and is not able to delay an event "sharp". We can do the same thing we did for spatialDistribution to fix that.

[Francesco Casella]

Replying to AnHeuermann:

Delay doesn't have a zero crossing function and is not able to delay an event "sharp". We can do the same thing we did for spatialDistribution to fix that.

That would be a good thing to do.

[Karim Adbdelhak]

I am not sure how the approach for concatination of zero crossings should work. There is some legacy code that collects them seperately, i am not sure if it is used at all.

@Francesco Does your approach hold for more than two conditions? So far i cannot find a reason it should not. I don't think it is worse if it contains discrete conditions, since it would behave the same if we are looking for a change in the discrete condition, but it would behave better on a change in the continous one.

We would need to map all logical operators to mathematical ones though. I don't know if it is that easy. If there is a way to do that there probably already exists a publication about it, we should search for one.

not a>0 -----> a<0
(a>0 and b>0) -----> (a*b)>0
(a>0 or b>0) -----> not(not a>0 and not b>0) 
             -----> not(a<0 and b<0)
             -----> not(-a>0 and -b>0)
             -----> not((-a)*(-b)>0)
             -----> (-a)*(-b)<0

Something like this. Most rules seem trivial but i think we want to represent everything with as few rules as possible to make it easy to maintain.

[Philip Hannebohm]

Replying to Karim.Abdelhak:

There you have it. If you compare

a > 2 and b > 3  --->  (a-2)*(b-3) > 0

with

a < 2 and b < 3  --->  (2-a)*(3-b) > 0  <=>  (a-2)*(b-3) > 0

they are indistinguishable! So this won't work in general.

I suggest, wee consider them separately or as a vector [a-2, b-3] and then check if it becomes all positive. In that case see which coordinate changed from negative to positive and use that for the bisection. If more than one changed, take the latest guess. And if it's a boolean just do interval halving as we do already.

[Francesco Casella]

Replying to Karim.Abdelhak:

I am not sure how the approach for concatination of zero crossings should work. There is some legacy code that collects them seperately, i am not sure if it is used at all.

I looked up Francois Cellier's book, but he doesn't get into this kind of detail. However, he only discusses continuous zcf's

[Francesco Casella]

Replying to phannebohm:

Replying to Karim.Abdelhak:

There you have it. If you compare

a > 2 and b > 3  --->  (a-2)*(b-3) > 0

with

a < 2 and b < 3  --->  (2-a)*(3-b) > 0  <=>  (a-2)*(b-3) > 0

they are indistinguishable! So this won't work in general.

@phannebohm, I'm not claiming that you should replace a < 2 and b < 3 with a single *condition* based on the product, which is of course ambiguous.

I just claim that you can define a single *zero-crossing function* to be passed to the solver, that contains the product, and can be used to detect the event times, when one of the original zcf's changes the sign. Whether or not that corresponds to an actual event, it needs to be checked with the original condition. If there is no event, e.g. a became < 2 and b was less than 3, then you can simply continue with the next integration step.

[Philip Hannebohm]

Replying to casella:

I just claim that you can define a single *zero-crossing function* to be passed to the solver, that contains the product, and can be used to detect the event times, when one of the original zcf's changes the sign. Whether or not that corresponds to an actual event, it needs to be checked with the original condition. If there is no event, e.g. a became < 2 and b was less than 3, then you can simply continue with the next integration step.

Okay, I was assuming the sign information of only the zcf is used by the solver to decide if it needs to find the zero or not. But then we need to do separate checks anyway afterwards. So I don't see the benefit.

Edit: You check one product-zcf and only if that changes sign you investigate further. Makes sense now!

[Francesco Casella]

Let's be very careful with scaling. If you combine p < p_max, W > Wmax, theta < theta_max, where p is a pressure, hence o.o.m. = 1e7 Pa, W is a power, hence o.o.m. = 1e9 W for a big power plant, and theta is a per-unit valve opening, hence o.o.m = 1 p.u., you may get numerically unreliable behaviour.

It is thus essential to use nominal attribute to rescale variables and zcf's. See our ​EOOLT paper on this topic.